{'fraction_non_alphanumeric': 0.044890114058825016, 'fraction_numerical': 0.0120887180395033, 'mean_word_length': 4.116056410919912, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 31, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Persistence and stability of generalized ribosome flow models with time-varying transition rates

In this paper the qualitative dynamical properties of so-called generalized ribosome flow models are studied. Ribosome flow models known from the literature are generalized by allowing an arbitrary directed network structure between the compartments and secondly, by assuming a general time-varying rate function describing the compartmental transitions. Persistence of the dynamics is shown using the chemical reaction network (CRN) representation of the system. We show the stability of different compartmental structures including strongly connected ones with an entropy-like logarithmic Lyapunov function. The L1 contractivity of solutions is also studied in the case of periodic reaction rates having the same period. It is also shown that different Lyapunov functions may be assigned to the same model depending on the factorization of the reaction rates.


Introduction
Compartmental models are used to describe and analyze the transport between different containers, called compartments in various natural and technological systems [1,2]. Compartments can be assigned to tissues or organs in pharmacokinetic models, mass containers in process systems, distinct disease states in epidemiological models, road sections in transportation systems or different habitats in ecological models. The modeled objects (molecules, people, vehicles, etc.) can move between compartments obeying the given constraints such as limits of directions, flow rates, or capacities. A fundamental feature of compartmental models is that each modeled object can be present in exactly one compartment at a given time. Naturally, compartmental models written in the original physical coordinates belong to the class of nonnegative systems for which the nonnegative orthant is invariant with respect to the dynamics [3,4]. This special property supports the dynamical analysis and control design in several ways. The controllability, observability, realizability and identifiability of mainly linear compartmental system are addressed in [5]. An excellent overview of the qualitative dynamical properties of general compartmental systems can be found in [6].
The dynamical modeling of the mRNA translation process has been in the focus of research since the second half of the 20th century (see, e.g. [7][8][9]). The first large scale analysis of gene translation through the so-called ribosome flow model (RFM) was presented in [10], where the applied second order nonnegative and nonlinear model based on the principle of Totally Asymmetric Exclusion [11] was able to capture the most important dynamical features of the translation process. In [12] the RFM was equipped with an appropriate input-output pair, and it was shown that after applying an affine positive output feedback, the system had a unique equilibrium point which is globally stable in the bounded operating domain. A circular RFM structure was analyzed in [13], where the authors proved using the theory of cooperative systems that the system has a continuum of equilibria, but each equilibrium is globally asymptotically stable within the equivalence classes of trajectories determined by the initial conditions. The stability of periodic solutions was also shown. In [14] a bounded pool of free ribosomes was added to the RFM generating a competition among the arbitrary number of mRNA molecules for ribosomes. This generates a special network structure for RFM subsystems, for which the uniqueness and stability of equilibria together with the properties of periodic solutions were also shown. Different compartment sizes of the RFM were assumed in [15], and it was shown that this modification does not change the favorable dynamical properties of the system. In [16], the ribosome flow model with Langmuir kinetics (RFMLK) is introduced, and a network structure is constructed with RFMLK subsystems connected through a pool. Among other results, it is shown that the trajectories of such a network always converge to a unique equilibrium.
Chemical reaction networks (CRNs) also called kinetic systems can be considered as universal descriptors of nonlinear dynamics, especially that of nonnegative systems [17]. Since the 1970's the theory of CRNs has been intensively studied, and there are several fundamental results on the relation between network structure/parametrization and dynamical properties [18]. The stability of mass-action type CRNs is most often analyzed using an entropy-like logarithmic Lyapunov function, originally called a "pseudo-Helmholtz function" in [19]. Probably the most well-known conjecture of chemical reaction network theory is the "Global attractor conjecture" according to which complex balanced kinetic systems are globally stable with respect to the nonnegative orthant with the logarithmic Lyapunov function [20]. This conjecture was proved for complex balanced reaction networks with a reaction graph of one component [21]. One of the most important results from the point of view of this paper is [22] studying zero deficiency networks, where the allowed kinetics is more general than mass action, the rate coefficients can be time-varying, and the logarithmic Lyapunov function is also generalized. The Lyapunov-function-based stability analysis of RFMs is mentioned as an important open problem in [23], which will be addressed in this paper using the CRN representation of the system.
It is interesting to mention that mathematical models which are equivalent to RFMs can also be obtained through a special finite volume spatial discretization of widely used flow models in PDE form [24]. These models also have a transparent representation in CRN form supporting further dynamical analysis. An arbitrary directed graph structure of such models with general time-invariant kinetics was considered in [25], where the existence and uniqueness of equlibria, persistence and contractivity (non-expansive property) of the solutions was shown using the theory of Petri nets, compartmental systems, and earlier results on RFMs. The stability of this model class with logarithmic Lyapunov functions was shown in [26], while a port-Hamiltonian description was given in [27].
Based on the above overview, the aim of this paper is to extend the results of [25] and [26] in the following respects: considering even more general kinetics with explicit timedependence, the qualitative analysis of periodic solutions, and finally, stability analysis with a family of different logarithmic Lyapunov functions.
The structure of the paper is the following. Section 2 contains the applied mathematical notations for compartmental models and kinetic systems. In Section 3 the kinetic representation of the studied model class is described, while new results on persistence and periodic behaviour in the time-varying case are proposed in Section 4. Stability analysis results with a family of non-unique logarithmic Lyapunov functions are described in Section 5, and finally, Section 6 summarizes the main results of the paper.

Notations and background
In this section, we describe the basic notations and building blocks of a compartmental system class and chemical reaction networks (CRNs). The notations and overview in this section are based on [25] and [27].

Compartmental models
Throughout the paper we consider systems containing a set of interconnected compartments and objects (such as ribosomes, particles, molecules, vehicles etc.) moving between them. We assume that the rate of transfer between compartments depends on the amount of objects in the source compartment as well as on the amount of free space in the target compartment. This naturally implies that each compartment has a well-defined finite capacity that limits the amount of modeled quantities that can be contained in the given compartment. We also allow explicit time dependence and in some cases dependence on the amount of objects and free space in other compartments.
For the formal definition, let us consider the set Q = {q 1 , q 2 , . . . , q m } of compartments and the set A ⊂ Q × Q of transitions, where (q i , q j ) ∈ A represents the transition from compartment q i into q j . Then, the directed graph D = (Q, A) is called the compartmental graph and it describes the structure of the compartmental model. The transitions are assumed to be immediate, thus loop edges are not allowed in the model since they do not introduce additional dynamical terms. Similarly, we do not allow parallel edges between two compartments in the same direction since they can be replaced by a single transition. We say that a (compartmental) graph is strongly connected if there exists a directed path between any two vertices in both directions, and we say that a graph is weakly reversible if it is a collection of isolated strongly connected subgraphs.
For each compartment q i we introduce the sets of donors and receptors, respectively, as that is, the set of donors of a given compartment are the compartments where an incoming transition originates from and the set of receptors are the compartments where an outgoing transition terminates in.

Chemical reaction networks (kinetic systems)
In this subsection we give a brief introduction of kinetic systems based on [18,19], where more details can be found. A chemical reaction network (CRN) contains a set of species Σ = {X 1 , X 2 , . . . , X N } and the corresponding species vector is given by X = [X 1 X 2 . . . X N ] T . The species of a CRN are transformed into each other through elementary reaction steps of the form where C j = y T j X and C j = y T j X are the source and product complexes, respectively, the vectors y j , y j ∈ N N 0 are stoichiometric coefficient vectors and functions K j : R N + × R + → R + are the rate functions with R + denoting the set of nonnegative real numbers. The matrix Y containing the stoichiometric coefficient vectors as columns is called the stoichiometric matrix. The subspace S ⊂ R N spanned by the so-called reaction vectors y j − y j is called the stoichiometric subspace of the CRN.
The CRN structure can be uniquely described by a directed graph as follows. For each complex we assign a vertex in the graph and and for each elementary reaction step of the form C j → C j we assign a directed edge between the corresponding vertices. We call the resulting graph the reaction graph of the CRN. The deficiency of the CRN is defined as δ = m − − s, where m is the number of distinct complexes, is the number of linkage classes (graph components) in the reaction graph and s is the dimension of the stoichiometric subspace.
Let x(t) ∈ R N + denote the state vector of the species as a function of time for t ≥ 0. Based on the above, the dynamics of the CRN is given bẏ We assume that a reaction can only take place if each species of the given reaction have nonzero concentration; that is, we assume that K j x(t), t = 0 whenever there exists k ∈ supp(y j ) such that x k (t) = 0, where we say that k ∈ supp(y j ) if [y j ] k > 0. This property ensures the invariance of the nonnegative orthant (or a part of it). We also presume standard regularity assumptions of the rate functions that guarantee local existence and uniqueness of solutions. Different results in this paper require different sets of such assumptions, thus for the sake of generality they will be specified later. Dynamics of the form of (2.2) is called persistent if no trajectory that starts in the positive orthant has an omega-limit point on the boundary of R N + . We note that for any v ∈ S ⊥ (where S denotes the stoichiometric subspace) we have that and thus x, v is constant. Since v ∈ S ⊥ was arbitrary we have that x(t) ∈ x(0) + S. This shows that the translates of S define invariant linear manifolds for the system. We further define for each p ∈ R n + a positive stoichiometric compatibility class S p = (p + S) ∩ R n + . A set of ODEs of the formẋ = f (x, t) is called kinetic if it can be written in the form (2.2) with appropriate rate functions and stoichiometric coefficient vectors.

Kinetic representation
In this section we construct a kinetic representation of the above compartmental system class. To do so, we assign a CRN that incorporates the compartmental structure. This allows the introduction of a system of ODEs of the form (2.2) describing the time evolution of the compartmental model. Some of the following steps are described in [25] or [27] in a time-invariant setting but here we recall and extend them for convenience.

Kinetic modelling of compartmental transitions
Let us consider a compartmental model D = (Q, A). Let the set of species be Σ = {N 1 , N 2 , . . . , N m } ∪ {S 1 , S 2 , . . . , S m } where N i and S i represent the number of particles and available spaces in compartment q i , respectively. To each transition (q i , q j ) ∈ A we assign a reaction of the form where K ij is the rate function of the transition. Such a reaction represents that during the transition from compartment q i to compartment q j the number of items decreases in q i and increases in q j , while the number of available spaces increases in q i and decreases in q j . Let n i and s i denote the continuous amount of particles and free space in q i , respectively.
Based on (2.2) the dynamics of the system is given bẏ where n and s denote the vectorized form of the variables n i and s i , respectively. It is easy to check that the model class in Eq. (3.1) contains ribosome flow models described in [23] or [15], and extends them in two ways: firstly, the reaction rate function K is not necessarily mass-action type and moreover, is time-varying, and secondly, the compartmental graph of the system can be arbitrary (i.e., there can be particle transition between any two compartments). Note, that we also allow the transition rates to depend on the amount of objects and free space in other compartments as well, perhaps describing inhibitory phenomena. Therefore, we call (3.1) a generalized time-varying ribosome flow model.
Clearly the reaction graph of the assigned CRN of a compartmental model is generally not strongly connected nor weakly reversible even if the compartmental graph is strongly connected. In fact, the reaction graph is weakly reversible if and only if each transition in the compartmental system is reversible. Even though the reaction graph, in some sense, loses the regularities of the compartmental graph, we can explicitly determine its deficiency from the compartmental topology and, as described in [25], CRNs of the form (3.1) exhibit persistence and stability properties in various senses in the time-invariant case.  Proof. For each transition between q i and q j we assign two complexes, namely N i +S j and S i + N j , regardless of the transitions' direction, so reversible reactions do not introduce additional complexes, and thus the number of stoichiometrically distinct complexes is M = 2|Ã|. A complex of the form N i + S j is only connected with the complex S i + N j , and thus we have = |Ã| linkage classes each consisting of exactly two complexes. To find the dimension of the stoichiometric subspace, denoted by S = dim S, observe that the reaction vector of a reaction of the form

Deficiency of CRNs realizing compartmental models
where e k ∈ R 2m denotes the kth unit vector. Again, since y i→j = −y j→i it suffices to consider the undirected graph |D|. Assume that y i→j is such that Then by (3.2) we have that for each non-zero term of the form c .→l y .→l the right-hand side also contains at least one non-zero term c l →. y l →. , including the terms c i→. y i→. and c .→j y .→j . This shows that the edges corresponding to the reaction vectors of the right-hand side form possibly multiple cycles in |D|. Without the loss of generality we may assume that this subgraph does not contain cycles isolated from (q i , q j ). We have to consider the following cases: 1. First, we assume that the right-hand side is a single chordless cycle and contains the transitions Taking the inner product of unit vectors e i , e l 1 , e l 2 , . . . , e lr , e j and c l k →l k+1 y l k →l k+1 + c lr→j y lr→j yields the system of linear equations: which clearly has one solution where each weight is equal to one.
2. If the right-hand side consists of multiple cycles, then repeatedly using the previous argument we can replace the arcs not containing (q i , q j ) with chords. Note, that if the reaction vector corresponding to the chord is already on the right-hand side, then we just have to modify its coefficient. This method decomposes the right-hand side and will leave us with one chordless cycle containing (q i , q j ), leading back to the previous case with exactly one solution. Repeating the arc substitutions we can see that each arc becomes a chordless cycle with the reintroduced edges and the arising systems of linear equations have exactly one solution.
The first case above shows that the dimension of the stiochiometric subspace reduces by one for each set of reaction vectors that correspond to edges forming a chordless cycle in |D| and the second case shows that is reduced by that exact amount. If σ denotes the number of chordless cycles inQ, then the deficiency of the reaction network can be computed as

Linear conservation laws
System (3.1) exhibits conservation in several senses. First of all, we have that thus the sum of modeled quantities and free spaces in the system is constant along the trajectories of (3.1); that is, the function H : is a first integral, where x 1 , x 2 , . . . , x m and x m+1 , x m+2 , . . . , x 2m correspond to the variables n 1 , n 2 , . . . , n m and s 1 , s 2 , . . . , s m , respectively. Our next observation is thatṅ i +ṡ i = 0 holds for each compartment, thus c i := n i + s i is the constant capacity of compartment q i . Let c (m) be a vector such that its ith coordinate is c i . Substituting s = c (m) − n we can rewrite (3.1) in a reduced state space aṡ or after an analogous substitution, aṡ As a consequence of the preceding observations, the functionH : is a first integral for (3.3), in which case each x i = n i (and similarly for (3.3) if each x i = s i ). This shows that while the state space of the decomposed systems isC : , for a given initial condition x(0)C the trajectories are contained in the (m − 1)-dimensional manifold (hyperplane) defined by For a generalized ribosome flow define c = n i=1 c i and for r ∈ [0, c] let L r ⊂C be the level set of H corresponding to r; that is, L r = a ∈C : H(a) = r .

Example 1.1
As a small example let us consider the compartmental model given by D = (Q, A), where The topology is shown in Figure 1. The corresponding CRN has the following species and reactions: It is easy to see that, indeed, the reaction graph is not weakly reversible and its deficiency is one. The dynamics of the model in the full state space is given by (3.1) aṡ which can be rewritten in the reduced state space based on (3.3) aṡ

Qualitative dynamical analysis
In this section we show that systems of the form (3.1) exhibit various interesting dynamical properties that can be characterized under different assumptions of the transition rate functions. First we will consider time-invariant systems to demonstrate the regularity of equilibria. Then we return to time-varying systems to generalize the results of [25].

Equilibria of time-invariant systems
In this subsection we assume that the K ij (n, s, t) rate functions are continuously differentiable and only depend on the variables n i and s j in a nondecreasing manner; that is, we assume that K ij (n, s, t) ≡ K ij (n i , s j ) for each i and j. Then the results [25,Propositions 5.5,5.6] show that a system of the form (3.3) is cooperative (the name also highlights the importance of the exclusion of inhibitory phenomena), is (strongly) monotone and each level set L r contains a unique globally (relative to its level set) asymptotically stable steady state. This implies that the steady states form a linearly ordered set. where n(0) ∈ L r is arbitrary and ρ t, n(0) denotes the solution at time t with ρ 0, n(0) = n(0). Clearly each e i is continuous and the monotonicity of the system also shows that each e i function is strictly increasing; that is, they are differentiable almost everywhere and their derivative are positive.

Example 1.2
Let us consider the triangular compartmental model in Figure 1. Its time evolution in the reduced state space is given in the form (3.3) aṡ and for the simulations we set capacities c 1 = 5, c 2 = 25, c 3 = 50. The rate functions in the different cases are assumed to have the form K ij (n i , c j −n j ) = k ij n i (c j −n j ) (corresponding to mass-action kinetics) or to be rational functions of the form for some l > 0 with k 12 = 100, k 23 = 40, k 31 = 60. Figure 2 shows the equilibrium curves for these rate functions with various l values. The topology is shown in Figure 3. The corresponding CRN has the following species and reactions: The dynamics of the system in the reduced state space is given bẏ Since the compartmental graph is not strongly connected the persistence and stability results of [25] are not applicable. However, empirical results show that the long-time behaviour of the system still exhibits some regularity, which can be divided into two cases base on the initial values of the system: and n 1 (t) and n 2 (t) will converge to the unique equilibrium on the level set . Note that since D is strongly connected, the results of [25] and the above investigation can be applied.
For the simulations we set c 1 = c 2 = c 3 = 100. The rate functions in the different cases are assumed to have form K ij (n i , c j − n j ) = k ij n i (c j − n j ) (corresponding to mass-action kinetics) or to be rational functions of the form for some l > 0 with k 23 = 15, k 32 = 25, k 31 = 35. Figure 4 shows the equilibrium curves for these rate functions with various l values. As described by the above cases we see that until the sum of the initial value exceed the capacity of the q 1 compartment the equilibrium lies on the n 1 axis. After that the equilibrium lies on the plane n 1 = c 1 ⊂ R 3 and since D is strongly connected we have that the coordinate functions of the equilibria e 2 (r) and e 3 (r), restricted to the set [c 1 , c], are continuous and strictly increasing. We note that while the system is not strongly connected it exhibits many similar qualitative properties as strongly connected models. For example, for initial values satisfying H n(0) > c 1 the system is Lyapunov stable as described in [27].  Based on the initial value and the exact compartmental structure the following phenomena can be observed: • Traps (and only traps) can become full, thus possibly creating new traps.
• Sources (and only sources) can become empty, thus possibly creating new sources.
• After a sufficient number of traps are filled and sources are emptied, the compartmental graph D is decomposed into isolated strongly connected components; that is, the resulting graph is weakly reversible, in which case the results of [25] can be applied.
While these observations are elementary and show that the system is stable, the equilibria are clearly non-unique with respect to the total mass of the network and in general it is not straightforward to predict from the initial value which components will fill and empty.

Persistence
In this subsection we consider time-varying generalized ribosome flows of the form (3.1) only under mild regularity assumptions described by the following theorem, which is based on the results of [28] but the statements are rephrased to be more aligned with our framework. For the definition of notions related to Petri nets (e.g. siphons) and their exact connection with CRNs we refer to [25,28]. (ii) There exists a positive linear conserved quantity c T x for the dynamics.
(iii) There are nonnegative, continuous functions K j (x), K j (x) such that (a) if x k >x k for each k ∈ supp(y j ), then K j (x) > K j (x) (and similarly for K j ) holds for each j = 1, 2, . . . , R, and (b) for each j = 1, 2, . . . , R, for all x ∈ R N + and for all t ≥ 0 we have To verify condition (i) we would, in general, need to enumerate all siphons of the CRN, which is well-known to be an NP-hard problem. However, in our recent paper [25] we explicitly characterized the siphons of a CRN assigned to a strongly connected compartmental models in the time-invariant case. However, one can observe that conditions (i) and (ii) of 4.2 are independent of the choice of transition rates and even independent from whether the system is time-invariant or not; that is, our results, formulated in the following theorem, hold for time-varying compartmental systems as well. It is not straightforward to determine exactly what types of reaction rates satisfy condition (iii). For the sake of specificity, we characterize a class of reaction rates of special interest which can be written in the following form where we assume that the transformations θ i , ν j ∈ C 1 (R) are nondecreasing, have θ i (0) = ν j (0) = 0 and satisfy 1 0 | log θ i (r)| dr < ∞ and 1 0 | log ν j (r)| dr < ∞ for each i, j = 1, 2, . . . , m. We also assume that the functions Ψ ij take the form Ψ ij (n, s) = α r (1) ,r (2)  where r (1) , r (2) ∈ N m and α r (1) ,r (2) ∈ R + . We further assume that for k ij (t) there exist k ij , k ij > 0 such that k ij ≤ k ij (t) ≤ k ij for all t ≥ 0. In this case we have which are clearly monotonous in the sense of Theorem 4.2, and thus condition (i) is satisfied and the system is persistent.
Remark 4.4. The above investigation and, in particular, condition (iii) of Theorem 4.2 shows that Lemmata 5.1, 5.2 and Remark 5.3 of [25] can be modified to the time-varying case; that is, for a system of the form (3.1) with strongly connected compartmental graph and reaction rates of the form (4.2), for each τ > 0 there exists holds for each i = 1, 2, . . . , m and t ≥ τ .
The denominator of (4.2) contains positive terms which can be interpreted as the inhibitory effect of other species, and the time-varying coefficient k ij (t) introduces the dependence of the system parameters on various factors such as temperature or the dynamical behaviour of other species that are not explicitly modelled as state variables. This class of rate functions contains many well-known examples, demonstrating the range and flexibility of reaction rates of the above form: 1. Setting each θ i (n i ) = n i and ν j (s j ) = s j and Ψ ij (n, s) = 0 we obtain the case of classical mass-action kinetics with time-varying rate coefficients: K ij (n, s, t) = k ij (t)n i s j .
2. Setting each θ i (n i ) = n i and ν j (s j ) = s j and Ψ ij (n, s) = l 2 − 1 + ln i + ls j + n i s j for some l > 0 yields K ij (n, s, t) = k ij (t) n i s j (l + n i )(l + s j ) corresponding to simple saturating kinetics described by the Monod equation.
3. The previous example can also be obtained by setting θ i (n i ) = n i l+n i and ν j (s j ) = s j l+s j and Ψ ij (n, s) = 0, showing that (4.2) is not unique. Notice however, that for fixed θ i , ν j transformations the function Ψ ij , and thus the fraction itself, is unique.
As a comparison let us consider the solutionñ(t) of the time-invariant system with the above nominal values. Figure 5a shows the phase portrait of the perturbed and the original systems starting from various initial conditions with H n(0) = 150. Figure 5b shows the time evolution of the state variables with n(0) = [5 45 100] T , where the state variables of the perturbed and the time-invariant system are depicted with blue lines and red lines, respectively. We can observe that since the time dependent terms are exponentially decaying and both systems evolve on the same linear manifold, the systems tend to the same equilibrium, as expected.

Periodic behaviour
Let us consider a generalized ribosome flow in the reduced state space of the form (3.3) with transition rates of the form (4.2) and assume that the transition functions are C 1 and periodic with the same period (but possibly different phase). Write (3.3) asṅ = F (t, n) and assume that the right-hand side satisfies the following monotonicity condition: F i (t, x) ≤ F i (t, y) for any two distinct points x, y ∈C such that x i = y i and x j ≤ y j for j = i. This condition is satisfied if, for example, the transition rates are such that Ψ ij ≡ 0; that is, if there are no inhibitory phenomena. Then the system phase locks (or entrains) with the periodic excitations.
Theorem 4.5. Consider a system of the form (3.3) satisfying the above monotonicity assumption, where each K ij (t) is periodic with a common period T . Then for each r ∈ [0, c] there exists a unique periodic function φ r : R + : →C with period T such that for all a ∈ L r we have that lim t→∞ ρ(t, a) − φ r (t) L 1 = 0.
Proof. The properties of the rate functions and the fact that ∇H is positive implies the result via [29], [30].
Remark 4.6. Since, in some sense, time-invariant systems can be seen as periodic, the stability result [25,Proposition 5.5] is a special case of the above theorem, where φ r is reduced to a single point of the manifold L r .

Stability analysis
In this section we show that generalized ribosome flows with reaction rate functions of the form (4.2) with piecewise locally Lipschitz k ij (t) coefficients satisfy a certain notion of robustness to the changes in the time-varying rate functions that can be traced back to the input-to-state stability of rate-controlled biochemical networks thoroughly investigated in [22]. The main difficulty in applying these results lies in the aforementioned fact that the CRN assigned to a compartmental model is generally not weakly reversible and its deficiency is generally not zero (see, Theorem 3.1) even if the compartmental topology is strongly connected. In order to circumvent this, we will perform a model reduction and rewrite (3.3) by factoring out certain terms. Let us first state the most important notions and results of [22].
Consider the system corresponding to a CRN with R reactionṡ where the nonnegative functions u ij are piecewise locally Lipschitz with a finite number of discontinuities and the stoichiometric coefficient vectors y i , y j are as described in 2.2.
In this section, however, we restrict the conditions on the transformation functions θ i : . Namely, we assume that (a) θ i is real analytic, Before continuing with the definitions, we consider the case when u(t) is a constant matrix A. We assume that A has nonnegative entries and is irreducible; that is, the underlying reaction graph is strongly connected. We denote the set of such A matrices as A. Then the equilibria ofẋ = f (x, A) can be divided into the sets of boundary equilibria and positive equilibria: Then, the result [22, Theorem 2.1] (and also [31,Theorem 2]) shows that if there are no boundary equilibria in any positive class, then each positive class contains a unique globally (relative to the positive class) asymptotically stable positive equilibrium. Denote the unique positive equilibrium in the same class as x 0 as x(x 0 , A) and notice that E A,+ = x(x 0 , A) x 0 ∈ R n + . Finally, denote Definition 5.1. We define the following function classes: (i) A function α : R + → R + is said to be of class K if it is continuous, strictly increasing and has α(0) = 0.
(ii) The subset of unbounded functions of class K are denoted by K ∞ .
(iii) A function β : R + × R + is said to be of class KL if β(., t) is of class K for all t ≥ 0 and β(r, .) is strictly decreasing to zero for all r > 0.
We consider nonnegative time-varying inputs such that at any time instant the reaction graph is strongly connected; that is, the input-value set U is a subset of A. Furthermore, let . 2 denote the spectral norm induced by the Euclidian norm and for u : u) is uniformly input-to-state stable (ISS) with inputvalue set U if for every compact set P ⊂ E and every compact set F ⊂ R n + containing P , there exist functions β = β P of class KL and φ = φ P of class K ∞ such that, for every x o ∈ P ∩ E u 0 ,+ for some u 0 ∈ U we have that holds for each u : R + → U input and every initial condition x 0 ∈ F ∩ S x 0 and for all t ≥ 0 such that x(s) ∈ F for s ∈ [0, t].
According to the above definition we say that a system is ISS if it is globally asymptotically system in the absence of external inputs and if its trajectories are bounded by an appropriate function of the input. In some sense this definition is intended to capture the idea of "bounded input bounded output" stability, since for bounded u input (u − u 0 to be more precise) the trajectories will remain in a ball and, in fact, approach the ball φ( u − u 0 U ) as t increases [32].
We assume that there exists a uniform lower bound on the parameters; that is, we consider input-value sets of the form We also recall that the input functions are piecewise locally Lipschitz in time with a finite number of discontinuities, thus we introduce Then the main Theorem of [22] states: Theorem 5.3. Consider the system (5) with and suppose that is is mass-conservative; that is, there exists v ∈ R n + such that v T f (x, u) = 0 for all x ∈ R n + and u ∈ A. Then the system with input maps u ∈ W is uniformly ISS with input-value set U ε .
The proof relies on the candidate ISS-Lyapunov function (for the definition of which and for the exact connection with ISS stability we refer to [22]) which, for mass-action systems, yields the classical entropy-like Lyapunov function wellknown from the theory of chemical reaction networks, see (5.4). We note that V (x, x) is uniquely determined by the θ i functions and does not depend explicitly on the reaction/compartmental structure or the time-varying u ij (t) functions; that is, it is universal in the sense of [33].
Remark 5.4. We note that the assumption that the compartmental graph (and thus the reaction graph of the factored model) is strongly connected is purely technical. For timeinvariant systems it simply ensures that the unique equilibrium on each level set of the first integral is positive (except for the trivial case of an empty network of course). In fact, in some cases the initial values of the network can ensure the positivity of the equilibrium even for not strongly connected systems (see Example 2 and [27] for more details), in which case the above Lyapunov function can be applied.

Factorization approach
Let us consider a generalized ribosome flow in the reduced state space of the form (3.3), in this case given bẏ Notice that we can naturally factor some terms of the transition rates into the timevarying coefficient as Then (5.1) can be rewritten aṡ This equation can be clearly embedded into the class of strongly connected systems of the form (5), since the reaction graph of (5.1) consists of species Σ = {N 1 , N 2 , . . . , N m }, has the m×m identity matrix as its stoichiometric matrix and for each transition (q i , q j ) ∈ A we assign a reaction of the form and thus the system of differential equations can be written aṡ where the elements ofÃ k are given by Note that the fractions 1+Ψ ij (n,c (m) −n) are differentiable (and thus Lipschitz) and each k ij (t) is piecewise locally Lipschitz, hence eachk ij (t) is piecewise locally Lipschitz. This shows that generalized ribosome flows can be embedded into the class of rate-controlled biochemical networks of [22] in a way that preserves the compartmental structure; that is, the reaction graph of (5.1) is topologically isomorph to the compartmental graph. In particular if the compartmental model is strongly connected, then the reaction graph of the reduced system (5.1) is strongly connected as well. Furthermore, combining the persistence of the system with Remark 4.4 we find thatÃ k ∈ W, and thus Theorem 5.3 ensures input-to-state stability.

Quasi-LTV factorization
A classical argument shows that the model reduction above can result in a Linear Time-Varying (LTV) system [6]. Consider an F (x) ∈ C k (R) nonnegative function such that F (0) = 0, where k ≥ 1. Then for the function F (rx) we have and thus and since F (0) = 0, we find that F (x) = xf (x). Note, that the calculation also shows that f ∈ C k−1 (R). Since θ i is real analytic we have that θ i (n i ) =θ i (n i )n i for someθ i real analytic function. Then (5.1) can be rewritten aṡ wherek Similarly as before, the reaction graph of (5.2) consists of species Σ = {N 1 , N 2 , . . . , N m }, has the m×m identity matrix as its stoichiometric matrix and for each transition (q i , q j ) ∈ A we assign a reaction of the form and thus the system of differential equations can be written aṡ where the elements ofÂ k are given by Again, eachk ij (t) is piecewise locally Lipschitz, thus for strongly connected compartmental models Theorem 5.3 ensures input-to-state stability via Remark 4.4.

Induced family of Lyapunov functions
The above investigation demonstrates that generalized ribosome flows can be embedded into rate-controlled biochemical networks in at least two different ways, where each embedding induces a different Lyapunov function of the form (5). Thus, in general, we may use at least two different Lyapunov functions governing the same dynamics. To characterize their exact connection, consider a factored system of the form (5.1) with its ISS-Lyapunov function V (n, n). The quasi-LTV representation of the system admits an ISS-Lyapunov function of the form While in general we are restricted to the above factorizations, in some special cases we may use a whole family of factorizations and corresponding Lyapunov functions. To illustrate this, consider an example when each θ i (r) = r a i (l+r) b i for some l > 0 and a i ∈ N, b i ∈ N 0 , a i ≥ b i (these properties ensure that the functions θ i are nondecreasing). Then, after the factorization described in Section 5.1, the Lyapunov function (5) becomes V (l,a,b) (n, n) = m i=1 (a i − b i )(n i − n i ) + a i n i log n i n i + b i (l + n i ) log l + n i l + n i .
We emphasize that (5) only depends on the θ i functions, in this case parametrized with the l, a i , b i values; that is, it is independent of the network structure and transition rate coefficients. We can also perform the factorization θ i (r) =θ i (r) râi (l+r)bi withâ i ∈ N,â i < a i , b i ∈ N 0 ,â i ≥b i yielding the Lyapunov function V (l,â,b) of the same form as in (5.4). This shows that the parameters a and b can be freely (apart from the constraints above) chosen in (5.4). We may also observe some interesting behaviour at the extrema of the parametersb and l, namely, that if we choose eachb i = 0 then the Lyapunov function in (5.4) is independent of l; that is, we have that where V LT V i is defined in (5.4).
function V (l,â,b) for various choices ofâ andb with l = 25 fixed. The second and third rows demonstrate the convergence characterized in (5.4); figures 7d-7f show V (l,â,b) for increasing l values 7g-7i shows m i=1â i V LT V i for the same increasing l values.

Example 3
Let us consider a compartmental system with m = 100 compartments in the reduced state space. We assume that the transition rate functions are corresponding to Hill kinetics (modified intentionally to have different powers in the numerator and the denominator) and are of the form with l = 350. We assume that the only nonzero coefficients are (n i − n i ) + 3n i log n i n i + n i log n 2 i + l n 2 i + l We can also factorize as θ i (r) =θ i (r) r 2 l+r 2 , when (5) becomes