Analysis of a reduced model of epithelial-mesenchymal fate determination in cancer metastasis as a singularly-perturbed monotone system
aa r X i v : . [ q - b i o . M N ] M a y Analysis of a reduced model of epithelialmesenchymalfate determination in cancer metastasis as asingularly-perturbed monotone system
M. Ali Al-Radhawi and Eduardo D. Sontag , Departments of Bioengineering and Electrical and Computer Engineering, , NortheasternUniversity, 360 Huntington Avenue, Boston MA 02115. Laboratory of Systems Pharmacology, Program in Therapeutic Science, Harvard Medical School,Boston MA 02115.Emails: { malirdwi,e.sontag } @northeastern.edu May 15, 2020
Abstract
Tumor metastasis is one of the main factors responsible for the high fatality rateof cancer. Metastasis can occur after malignant cells transition from the epithelialphenotype to the mesenchymal phenotype. This transformation allows cells to migratevia the circulatory system and subsequently settle in distant organs after undergoingthe reverse transition from the mesenchymal to the epithelial phenotypes. The coregene regulatory network controlling these transitions consists of a system made upof coupled SNAIL/miRNA-34 and ZEB1/miRNA-200 subsystems. In this work, weformulate a mathematical model of the core regulatory motif and analyze its long-termbehavior. We start by developing a detailed reaction network with 24 state variables.Assuming fast promoter and mRNA kinetics, we then show how to reduce our modelto a monotone four-dimensional system. For the reduced system, monotone dynamicalsystems theory can be used to prove generic convergence to the set of equilibria forall bounded trajectories. The theory does not apply to the full model, which is notmonotone, but we briefly discuss results for singularly-perturbed monotone systemsthat provide a tool to extend convergence results from reduced to full systems, underappropriate time separation assumptions.
Realistic dynamical models of physical systems are often very complex and high-dimensional,and this is especially so in molecular cell biology. Effective analysis often requires simplifi-cations through model reduction techniques. A traditional approach to model reduction inbiology is to take advantage of time-scale separation. Indeed, most interesting processes inbiology are made up of subsystems which operate at different time scales, thus allowing fast1ubprocesses to be “averaged out” at the observational time scale. The rigorous mathemat-ical analysis of time scale separation, and in particular the mathematical field of singularperturbations , owe much to the pioneering work of Tikhonov [1]. Singular perturbationtheory now plays a key role in both science and engineering [2].There are many examples of the ubiquitousness of time-scale separation in molecularbiology. An early example is the study of enzymatic reactions, where the derivation of
Michaelis-Menten kinetics in 1913 [3] is still widely used [4]. Another example is providedby Gene Regulation Networks (GRNs) which naturally have multiple levels of time scales:external stimuli change Transcription Factor (TF) activities in milliseconds, promoter kinet-ics equilibrate in seconds, transcription and translation take minutes, and protein kineticsare in the order of tens of minutes to hours [5].In this paper we study a GRN that determines cell-fate in the metastasis of canceroustumors. This network regulates a transition between two cell types: epithelial cells , whichline the external and internal surfaces of many organs, and mesenchymal stem cells , whichare multipotent connective tissue cells that can differentiate into other type of cells such asmuscles, bone, etc. Bidirectional transitions between these two cell types can happen, andthey are referred to as “Epithelial to Mesenchymal” transitions (EMT) and “Mesenchymal toEpithelial” transitions (MET). During the EMT process, a cell loses adhesion to neighbouringcells, and becomes more invasive and migratory. It is worth noting that both EMTs andMETs are part of normal developmental processes such as embryogenesis and tissue healing.Nevertheless, they are of one of main mechanisms of tumor metastasis. After undergoingEMTs, cancerous cells travel through the blood as Circulating Tumor Cells (CTCs). TheseCTCs settle in other organs by undergoing METs, and they subsequently multiply, thusgiving rise to metastatic tumors [6, 7]. Because of their key role in cancer, the identificationof the GRNs enabling EMTs/METs has been a focus of a large research effort. Transcriptionfactors (TFs) such as SNAIL, SLUG, TWIST, and ZEB1 have been studied in great detail[8, 9]. Cellular signals such as p54, Notch, EGF, Wnt, HIF-1 α , and others can induceEMTs/METs [6]. These signals act on a core four-component network that involves theupstream SNAIL/miR-34 circuit and the downstream ZEB1/miR-200 circuit [10, 11]. Theconceptual organization for such a circuit is depicted in Figure 1. Each pair reproducesthe standard toggle switch architecture, i.e., mutual inhibition. However, this design differsby having a mixed inhibition mechanism: the TFs ZEB1 and SNAIL inhibit the miRNAs( µ , µ ) at the transcriptional level, while µ , µ inhibit SNAIL and ZEB1, respectively,at the translational level. Therefore, such circuits are known as chimeric circuits [10].One of the contributions of this work is to translate the conceptual diagram in Figure 1into a precise mathematical model. Our model, while not completely novel, clarifies andprovides explicit details of several features of models found in the literature. More impor-tantly, we analyze the long-term behavior of a reduced model obtained under a natural andbiologically realistic time scale separation assumption. We prove a theorem guaranteeing“almost-sure” convergence of trajectories to steady states. This paper also serves to moti-vate a powerful but relatively unknown theorem on singularly perturbed monotone systems.A central and well-known result for monotone systems is Hirsch’s Generic Convergence The-orem [12, 13, 14, 15, 16], which guarantees that almost every bounded solution of a stronglymonotone system converges to the set of steady states. This theory has been widely appliedto biochemical systems [17, 18, 19]. However, many biological models are not monotone. For2igure 1: The core EMT/MET network [10], which involves proteins, mRNAs, microRNAs(miRNAs) and the underlying genes which are not explicitly depicted. The mRNAs ofSNAIL and ZEB1 are denoted by m s , m z , respectively. An arrow of the form “ → ” denotesactivation, while “ ⊣ ” denotes inhibition. The detailed reaction network model is presentedin § S and m S (seeFigure 1). Nevertheless, we can take advantage of the fact that transcription happens ata faster time-scale than translation. Moreover, despite the fact that miRNAs and mRNAsbelong to the same class of biochemical molecules, comparison of their half-lives reveals thatmRNAs have much shorter half-lives than miRNAs and proteins [20, 21]. Hence, the neg-ative loop in the core EMT/MET network is a “fast” loop. Intuitively, negative loops thatact at a comparatively fast time scale should not affect the main characteristics of monotonebehavior. This has been rigorously studied in [22] (see the Ph.D. thesis [23] for more details),using tools from geometric singular perturbation theory.The paper is organized as follows. In Section II we model the core EMT/MET network indetail, and derive a reduced model. Section III reviews several basic definitions and theoremsabout monotone systems, and in Section IV we analyze the reduced model theoretically andwe review results that are can be used to extend convergence from reduced to full systems. A mathematical model for the EMT/MET system has been presented in [24, 10] based onseveral assumptions that include detailed balance. Here, we drop such assumptions, anddevelop a more “mechanistic” model via the framework of Chemical Reaction Networks3CRN). We review the notation briefly [25, ? ]. A CRN is a set of species S = { Z , .., Z n } and a set of reactions R = { R , ..., R ν } . A reaction R j can be written as: P ni =1 α ij Z i → P ni =1 β ij Z j .A stoichiometry matrix Γ ∈ R n × ν is defined elementwise as [Γ] ij = β ij − α ij . The reactionsare associated with a rate function R : R n ≥ → R ν ≥ . We assume that R takes the form ofMass-Action kinetics: R j ( z ) = Q ni =1 k j z α ij j , where k j is the kinetic constant. Let z ( t ) ∈ R n ≥ be the vector of species concentrations at time t . The associated ODE can be written as:˙ z = Γ R ( z ).We model GRNs as CRNs via the central dogma of molecular biology (see [27] for adetailed framework). Each gene is associated with promoter states, an mRNA state, and aprotein state. As mentioned in the introduction, we can assume that the promoter kineticsand the mRNA kinetics are fast. In this section, we will derive the dynamics of the fastand slow systems. We will show that the model can be reduced from
24 dimensions to fourdimensions . The slow states are S, Z, µ , µ . We will start with promoter kinetics. For each gene j , we denote the promoter states by species of the form D ji . The super-script denotes the gene, while the subscript denotes the occupancy of the binding sites. Let D , D , D S , D Z be the unbound promoters for µ , µ , SNAIL, and ZEB1, respectively.For instance, D s denotes S binding to the promoter of µ , while D sz denotes both S, Z binding to the promoter of µ .Hence, the CRN describing promoter dynamics for the EMT/MET network (Figure 1)can be written as:2 S + D Z α s −− ⇀↽ −− α − s D Zs , Z + D Z α z −− ⇀↽ −− α − z D Zz , S + D Zz α s −− ⇀↽ −− α − s D Zsz , Z + D Zs α z −− ⇀↽ −− α − z D Zsz , S + D α s −− ⇀↽ −− α − s D s , Z + D α z −− ⇀↽ −− α − z D z , S + D z α s −− ⇀↽ −− α − s D sz , Z + D s α z −− ⇀↽ −− α − z D sz , S + D α s / −− ⇀↽ −− α − s D s , Z + D α z / −− ⇀↽ −− α − z D z , S + D z α s / −− ⇀↽ −− α − s D sz , Z + D s α z / −− ⇀↽ −− α − z D sz , S + D S α s / −− ⇀↽ −− α − s D Ss , The concentration of a species Y will be denoted, if no other notation is used, by [ Y ]. Thereaction structure implies that we have the conservation law P i [ D ji ]( t ) = E j for each gene,where E j , j ∈ { S, Z, , } is the total concentration available in the medium. Hence, E j stays constant during the course of the reaction. Note that each TF is assumed to bind toits promoter as a dimer. The same analysis can be repeated if the TF binds as an n -mer forsome integer n ≥ active promoter states. Since both S, Z repress µ , µ , then theactive states are the unbound promoter states D , D , respectively. The promoter D S is4he active state for SNAIL since it is a self-repressing gene. Finally, D Zsz is the active statefor ZEB1 because it is activated by both SNAIL and self-binding.We will consider ZEB1 as an example. Let d z ( t ) := [[ D Z ] , [ D Zs ] , [ D Zz ] , [ D Zsz ]] T ( t ). Usingthe reactions above, we can write the ODE for the promoters of Z to get:˙ d z = P z ( s, z ) d z := − α s s − α z z α − s α − z α s s − α − s − z α z α − z α z z − α − z − s α s α − s z α z s α s − α − s − α − z d z . (1)Note that s, z are slow variables. Setting the derivatives to zero and substituting the con-servation law we get: [ D Zsz ] qss = E Z s z ( s + A s )( z + A z ) , where A s := α − s /α s , A z := α − z /α z . Similarly, we can define the matrices P s ( s, z ) , P ( s, z ) , P ( s, z ).Hence, we get:[ D ] qss = E A s A z ( s + A s )( z + A z ) , [ D S ] qss = E S A s ( s + A s ) , [ D ] qss = E A s A z ( s + A s )( z + A z ) . Translational inhibition by miRNAs is achieved by an RNA-induced Silencing Complex(RISC). RISC binds to the target mRNA and degrades it via the Argonaute protein [28].Here, we assume a simple model in which miRNA binds to the target mRNA to inhibit trans-lation and activate degradation. In general, multiple miRNAs can bind to a single mRNA.Here we assume that each mRNA can have up to two binding sites. For each additionalmiRNA binding, translation is inhibited and degradation is accelerated.The CRN for the transcription of SNAIL and ZEB1 can be written as follows: D S γ s −→ D S + M s , D Zsz γ z −→ D Z + M z , M s β s −→ ∅ , M z β z −→ ∅ , where M s , M z are the species denoting the mRNAs of SNAIL and ZEB1, respectively. TheCRN for the miRNA-mRNA reaction (with two binding sites) can be written as follows: µ + M S c s −− ⇀↽ −− c − s M Sµ β s −−→ µ , µ + M Sµ c s −− ⇀↽ −− c − s M Sµ β s −−→ µ ,µ + M Z c z −− ⇀↽ −− c − z M Zµ β z −−→ µ , µ + M Zµ c z −− ⇀↽ −− c − z M Zµ β z −−→ µ , Since miRNAs actively degrade mRNAs, we assume that: β s < β s < β s , β z < β z < β z . (2)Remembering that µ , and z are the slow variables, we can write the ODE for ZEB1translational dynamics as follows (with m z ( t ) := [[ M z ] , [ M zµ ] , [ M zµ ]] T ( t )):˙ m z = Q z ( µ ) m z + b z D Zsz (3):= − c z µ − β z c − z c z µ − β z − c z µ − c − z c − z c z µ − c − z − β z m z + γ z D Zsz µ , the quasi-steady state can be solved easily by matrixinversion. Substituting the QSS for [ D Z ] we get: m z,qss = γ z E Z s z ( s + A s )( z + A z ) β z c z µ + e z e z µ c z µ β z c z µ + e z c z µ + e z β z , (4)where e z := β z β z + ( β z + β z ) c − z + c − z , e z := c z ( β z + c − z ) , e z := β z β z + β z β z + β z c − z .The dynamics of m s can be derived similarly and we get: m s,qss = γ s E S A s ( s + A s ) β s c s µ + e s e s µ c s µ β s c s µ + e s c s µ + e s β s , (5)where e s := β s β s + ( β s + β s ) c − s + c − s , e s := c s ( β s + c − s ) , e s := β s β s + β s β s + β s c − s . The complete reaction model requires modeling the production of the proteins and themiRNAs. The CRN can be written as follows: D εβ −−→ D + µ , µ εβ − −−−→ ∅ ,D εβ −−−→ D + µ , µ εβ − −−−−→ ∅ ,M s εk s −−→ M s + S, M sµ εk s −−→ M sµ + S, M sµ εk s −−→ M sµ + S, S εk − s −−→ ∅ ,M z εk z −−→ M z + Z, M zµ εk z −−→ M zµ + Z, M zµ εk z −−→ M zµ + Z, Z εk − z −−→ ∅ . The kinetic constants are multiplied by ε to emphasize that the corresponding reactionsare slow. Note that since miRNA inhibits translation, the mRNA-miRNA complexes havelower translation rates than raw mRNAs. Hence, we have: k s > k s > k s , k z > k z > k z . (6)The model above assumes that Z is produced only when both Z , and S are bound toZEB1’s promoter. But this is not realistic, since constitutive transcription is always presentat low levels. Hence, we model this by a reaction ∅ εδ z −→ Z with δ z small. This will help usalso show generic convergence for the reduced system in § µ = β [ D ] − β µ = E A z A s β ( s + A s )( z + A z ) − β µ ˙ s = [ k s k s k s ] m s − k − s s = [ k s k s k s ] m s,qss − k − s s (7)˙ µ = β [ D ] − β µ = E β A z A s ( s + A s )( z + A z ) − β µ ˙ z = δ z + [ k z k z k z ] m z − k − z z = δ z + [ k z k z k z ] m z,qss − k − z z. where m z,qss , m s,qss are given by (4),(5), respectively.6 Monotone systems and singular perturbations
In this section, we review basic definitions and results regarding monotone systems. We baseour discussion on [13, 16, 23].A nonempty, closed set C ⊂ R N is said to be a cone if C + C ⊂ C , αC ⊂ C for all α > C ∩ ( − C ) = { } . For each cone C , a partial order on R N can be associated. For any x, y ∈ R N , we define: x ≥ y ⇔ x − y ∈ Cx > y ⇔ x − y ∈ C, x = y. When C ◦ is not empty, we can define x ≫ y ⇔ x − y ∈ C ◦ .In this paper, we only consider cones which are orthants of R N . In order to identify thevarious orthants, let σ ∈ {± } N and let C σ = { z ∈ R n | σ i z i ≥ , ≤ i ≤ n } be the corresponding orthant cone. Let (cid:22) σ be the corresponding partial order.A set W ⊆ R N is said to be p -convex, if W contains the line joining x and y whenever x (cid:22) y , x, y ∈ W . Hence, equipped with a partial ordering on a p -convex and open set W westudy the ordinary differential equation: dzdt = F ( z ) , (8)where F : W → R N is a C vector field. We assume that the system is forward invariantwith respect to W . In our application in this paper we will use W = R N + , the open positiveorthant.We are interested in a special class of equations which preserve the partial order alongall trajectories. Definition 1.
The flow φ t of (8) is said to have positive derivatives on a set W ⊆ R N , if (cid:2) ∂∂z φ t ( z ) (cid:3) x ∈ C ◦ for all x ∈ C \ { } , z ∈ W , and t > . Definition 2.
The system (8) is called monotone (resp. strongly monotone) on a set W ⊆ R N if for all t > and all z , z ∈ W , z ≥ z ⇒ φ t ( z ) ≥ φ t ( z ) ( resp. φ t ( z ) ≫ φ t ( z ) when z = z ) . Establishing that a flow has positive derivatives can be performed by verifying the irre-ducibility of the Jacobian as the following theorem states:
Proposition 1.
Assume that (8) is monotone with respect to an orthant cone C . If ∂F∂z ( z ) isirreducible for all z ∈ W then the flow φ t of (8) has positive derivatives on the set W ⊆ R N . The following result can be interpreted as saying that having positive derivatives is an”infinitesimal” version of strong monotonicity.7 roposition 2.
Let W ⊂ R N be p -convex and open. If the flow φ t has positive derivativesin W , then the associated ODE is strongly monotone on W . Hence, we have the following corollary:
Corollary 1.
Let (8) and W be given as above. Assume that (8) is monotone with respectto an orthant cone C . If ∂F∂z ( z ) is irreducible for all z ∈ W , then the system (8) is stronglymonotone with respect to C on W . This is a rephrasing of the well-known Hirsch Generic Convergence Theorem:
Theorem 1.
Assume that (8) is strongly monotone on a p -convex open set W ⊆ R N . Let W c ⊆ W be defined as the subset of points whose forward orbit has compact closure in W . Ifthe set of equilibria is totally disconnected, then the forward trajectory starting from almostevery point in W c converges to an equilibrium. By Corollary 1, it is sufficient to check monotonicity and irreducibility to establish genericconvergence.
Monotonicity with respect to orthant cones can be characterized via
Kamke’s conditions . Wereview the relevant material from [16, 18]. Let σ ∈ {± } N and let C σ be the correspondingorthant as defined above. Let Σ = diag( σ ), i.e., a diagonal matrix with σ as the diagonal.Then, the following holds: Theorem 2. (Kamke’s conditions) Let (8) be given and let φ t be the associate flow. Let J = ∂F∂z be the corresponding Jacobian. Let W ⊆ R N be a p -convex open set. Assume thatthere exists σ ∈ {± } N such that Σ J Σ is Metzler on W (i.e., all non-diagonal entries arenon-negative). Then the corresponding flow is monotone on W with respect to the partialorder (cid:22) σ . The last proposition gives a useful characterization for monotonicity with respect toorthant cones; however, it requires checking 2 N possible sign combinations. Alternatively, asimple graphical criteria can be stated for a graph derived from the Jacobian. Informally, itstates that the system is orthant monotone if every loop has a net positive sign . We state itmore formally next.We say that a Jacobian J is sign-stable on a set W if for each i = j , sgn( J ij ) is constanton W . Define a signed directed graph G with vertices { , ..., n } . There is an edge connectingvertices i to j if the partial derivative J ij does not vanish identically. The sign of the edge isequal to sign of J ij . A loop is any sequence of edges (without regard to direction) that doesnot traverse a vertex twice and it starts and ends with the same vertex. The sign of a loopis the product of the signs of the constituent edges. It is worth noting that diagonal entries,i.e., “self-loops”, have no effect on the validity of the results. Theorem 3. (Positive Loop Property) Let (8) be given. Let J = ∂F∂z be the correspondingJacobian. Assume J is sign-stable. Let W ⊆ R N be a p -convex open set. Let G be the graphdefined above. If every loop has a positive sign, then there exists σ ∈ {± } N such that Σ J Σ is Metzler on W (i.e, all non-diagonal entries are non-negative). .3 Singular perturbation model reduction The process of mathematical modeling of physical processes involves usually reduction of thedynamics of fast variables. For instance, self-loops, i.e., terms corresponding to the diagonalentries of the Jacobian J , are often approximations of fast dynamics. Hence, it is intuitive toexpect that the theorems in the previous section hold for sufficiently fast negative self-loops.Hence, we consider a system in the singularly perturbed form: dxdt = f ( x, y, ε ) (9) ε dydt = g ( x, y, ε ) , where f : R n + × R m + × [0 , ¯ ε ] → R n + , g : R n + × R m + × [0 , ¯ ε ] → R m + are smooth bounded functionsand ¯ ε > ε < ¯ ε . These assumptions are automatically satisfied in our case since Mass-Actionkinetics give rise to polynomial systems.For 0 < ε ≪
1, the dynamics of x are much slower than y . If ε = 0, we can change thetime scale to τ = t/ε , and study the equivalent form: dxdτ = εf ( x, y, ε ) (10) dydτ = g ( x, y, ε ) . Model reduction via singular perturbations requires solving the fast system at a quasi-steady state, hence we assume that there exists a smooth bounded function m : R + m → R + n such that g ( x, m ( x ) ,
0) = 0 for all x ∈ R n + . From the previous section it can be seen that m exists and is unique as we have derived all the QSS expressions uniquely.For simplicity, let z = y − m ( x ). Hence, the fast system (10) can be written as: dxdτ = εf ( x, z, ε ) (11) dzdτ = g ( x, z, ε ) , where f ( x, z, ε ) = f ( x, z + m ( x ) , ε ) ,g ( x, z, ε ) = g ( x, z + m ( x ) , ε ) − ε [ ∂∂x m ( x )] f ( x, z, ε ) . When ε = 0, the system (11) becomes dzdτ = g ( x, z, , x ( τ ) ≡ x ∈ W x . (12)Verifying the stability of the fast system is a standard and an intuitive requirement forsingular perturbation methods as in the original Tikhonov theorem [1]. Hence we need tocheck that: 9. the equilibrium z = 0 of (12) is globally asymptotically stable on { z | z + m ( x ) ∈ R m + } for all x ∈ R + n .2. all eigenvalues of the Jacobian ∂∂y g ( x, m ( x ) ,
0) have negative real parts for every x ∈ R + n .We study this in the next section. In our analysis in the section 2 we have considered fast and slow reactions separately. We willuse now the notation from §
3. Let x ( t ) = [ µ , s, µ , z ] T ( t ) be the state of the slow system,and let y ( t ) = [ d Ts , d Tz , d T , d T , m Ts , m Tz ]( t ) be the state of the fast system. Hence, the fullsystem can be written in the form (9). We denote the reduced system (7) as ˙ x = G ( x ). In this subsection, we want to study global asymptotic stability of the fast system, and toshow that its Jacobian is Hurwitz.For a fixed slow variable x , the fast dynamics of ZEB1, SNAIL, miRNA-34 , miRNA-200are decoupled from each other as can be seen from the CRNs introduced before. Furthermore,all the systems are linear. Let us analyze ZEB1 as an example. From (1),(3) we get (cid:20) ˙ d z ˙ m z (cid:21) = (cid:20) P z ( s, z ) 0 B z Q z ( s, z ) (cid:21) (cid:20) d z m z (cid:21) , where B z := [ O × , b z ] and O × ∈ R × is a zero matrix. The dynamics of d z are decoupledfrom m z . Since P z is an irreducible Metzler matrix with principal eigenvalue of 0, Perron-Frobenius Theorem [29] implies that all other eigenvalues are strictly negative. In fact, theycan be computed as {− α s s − α − s , − α z z − α − z , − α s s − α − s − α z z − α − z } . Since thetotal number of promoters is conserved, we can reduce the dimension of the ODE by one. Itfollows that the reduced promoter dynamics are globally asymptotically stable. The sameargument applies to the promoters of Z, µ , µ .We turn to the dynamics of m z . Due to the block triangular structure we only need tostudy the eigenvalues of Q z ( s, z ). This matrix is also Metzler and the principal eigenvalue canbe shown to be negative. Hence, it is Hurwitz. This can be proven alternatively with a linearLyapunov function V ( m z ) = T m z . The time-derivative is ˙ V ( m z ) = − [ β z β z β z ] m z < m z = 0. The same argument applies to the translational dynamic of SNAIL. Hence,the fast system is globally asymptotically stable and its Jacobian is Hurwitz. We study the reduced system (7). We will utilize Kamke’s conditions as in Theorem 2.Hence, we compute the Jacobian J for the slow system (7):10igure 2: Graph of the reduced circuit. In the terminology of Theorem 3, “ → ” has a + sign,and “ ⊣ ” has a − sign. J = − β − E A s A z β s ( s + A s ) ( z + A z ) − E A s A z β z ( s + A s ) ( z + A z ) J J − E A s A z β s ( s + A s ) ( z + A z ) − β − E A s A z β z ( s + A s ) ( z + A z ) J J J , where J = − h ( µ ) E s A s γ s ( s + A s )( β s c s µ + e s c s µ + e s β s ) , where h ( µ ) := c s ( β s + c − s ) ( β s + c − s )( β s k s − β s k s )+2 c s µ ( k s β s − β s k s )( c − s + β s )( c − s + β s ) + c s µ ( β s ( β s k s − β s k s ) + ( β s + c − s )( β s k s − β s k s )). Note that h ( µ ) > J <
0. Similarly, we can also calculate J toshow that it is negative whenever (2),(6) are satisfied.Finally, using (4) we compute J as: J = 2 E Z γ z A s sz ( s + A s ) ( z + A z ) k z ( β z c z µ + e z ) + k z ( e z µ ) + k z c z µ β z c z µ + e z c z µ + e z β z > . Remark 1.
Conditions (2) , (6) are stronger than we need. In fact, it can be seen from theexpression of h above that it is sufficient to have the protein production ratio of the rawmRNA greater than the corresponding one for the miRNA-RNA complex. The signs of J , J do not play a role since they correspond to self-loops. Therefore, wehave the following sign pattern for the Jacobian under the conditions (2),(6): ∗ − −− ∗ − ∗ − − ∗ , σ ∈ {± } N such that Σ J Σ is Metzler.Therefore, by Theorem 2 we get that the reduced system (7) is monotone on W = R .Monotonicity holds with respect to the orthant cone specified by σ = [ − , , − , T . (see § open positive orthant R n + . Hence,by Corollary 1 the reduced system is strongly monotone on R n + .Fix an initial condition x ∈ R n + , and let x ( t ) := ϕ ( t ; x ) be the corresponding solutionof the slow system. In order to infer generic convergence to the set of equilibria, we needthe strong order preserving property to hold on the ω -limit set of the solution. Hence, weneed to show that ω ( ϕ ( t ; x )) ∩ ∂ R n + = ∅ . In other words, we need to show that the slowsystem is persistent . This can be shown as follows. Recall that the slow system (7) is givenas ˙ x = G ( x ). It can be verified from (7) that ∀ i, G i ( z ) > z i = 0. Hence, thereexists η > x i ( t ) > x i ( t ) ∈ [0 , η ). This implies that the boundaryhas no equilibria and is repelling (i.e., the vector field is pointing away from the boundary).We proceed by seeking contradiction. W.l.o.g, assume that ∃ x ∗ ∈ ω ( ϕ ( t ; x )) ∩ ∂ R n + with x ∗ i = 0 (the i -coordinate). Hence, there exists a sequence { t k } ∞ k =1 such that x i ( t k ) > x i ( t k ) → x ∗ i = 0. Therefore, for each right neighborhood N of 0 there exists t ∗ for which x i ( t ∗ ) ∈ N and ˙ x i ( t ∗ ) <
0, which is a contradiction. Hence, the reduced system is persistent.Therefore, we state the following theorem.
Theorem 4.
Consider the reduced system (7) . Let W c ⊂ R be the subset of points whoseforward orbit has a compact closure in R . Then, the forward trajectory starting from almostevery point in W c converges to an equilibrium. An interesting general question for singularly perturbed systems is as follows. Assumingthat the slow system is strongly monotone, does the full system obey generic convergenceproperties?
In other words, suppose that the flow ψ t of the slow system (set ε = 0 in (9)): dxdt = f ( x, m ( x ) ,
0) (13)has strong monotonicity properties that guarantee almost-global convergence, meaning con-vergence to equilibria for all initial states except for those states in a set of measure zero(or, in a topological formulation, a nowhere dense set), but that the complete system isnot monotone (so that no such theorem can be applied to it). Still, one may expect thatthe almost-convergence result can be lifted to the full system, at least for small ε >
0. Anobstruction to this argument is that there is no a priori reason for the exceptional set tohave a pre-image which has zero measure (or is nowhere dense). Notheless, a positive resultalong these lines was developed in [22] via the use of geometric invariant manifold theory toexamine the fibration structure and utilize an asymptotic phase property [30, 31, 32]. Figure3 illustrates the idea. The ODE restricted to the invariant manifold M ε can be seen as aregular perturbation of the slow ( ε =0) ODE. In [13], it has been noted that a C regular12igure 3: Illustration of the manifolds M and M ε . The figure shows two key properties of M ε . First, M ε is close to M . Second, the trajectories on M ε converge to steady-states ifthose on M do.perturbation of a flow with positive derivatives inherits the generic convergence properties.So, solutions in the manifold will generally be well-behaved, and asymptotic phase meansthat trajectories close to M ε can be approximated by solutions in M ε , and hence they alsoapproach the set of equilibria if trajectories on M ε do. We refer the reader to [22] for atechnical discussion. In principle, this approach could be applied to systems such as ours, toconclude almost global convergence, under mild technical conditions on domains of validityfor equation. We omit the details of the application here. We have conducted a model reduction analysis for the EMT/MET system. Bounded tra-jectories of the reduced EMT/MET system generically converge to steady states, assumingsufficiently fast promoter and mRNA kinetics. Therefore, such a model cannot admit oscil-lations nor chaotic behavior. There are many further directions for research, which are beingpursued by the authors, including the generalization to the case of miRNA binding to morethan two sites on the mRNA molecule.
Acknowledgements
We thank the reviewers of this manuscript for the meticulous reading and the constructiveremarks. This research was supported by NSF grants 1716623 and 1849588.
References [1] A. N. Tikhonov. Systems of differential equations containing small parameters in thederivatives.
Matematicheskii sbornik , 73(3):575–586, 1952.[2] P. Kokotovic, H. K. Khalil, and J. O’Reilly.
Singular perturbation methods in control:analysis and design . SIAM, 1999. 133] L. Michaelis and M. L. Menten. Die kinetik der invertinwirkung.
Biochem Z. , 49:333–369, 1913.[4] J. Gunawardena. Time-scale separation–Michaelis and Menten’s old idea, still bearingfruit.
The FEBS journal , 281(2):473–488, 2014.[5] U. Alon.
An introduction to systems biology: design principles of biological circuits .Chapman and Hall/CRC, 2006.[6] J. P. Thiery, H. Acloque, Ruby Y. J. Huang, and M. A. Nieto. Epithelial-mesenchymaltransitions in development and disease. cell , 139(5):871–890, 2009.[7] A. W. Lambert, D. R. Pattabiraman, and R. A. Weinberg. Emerging biological princi-ples of metastasis.
Cell , 168(4):670–691, 2017.[8] B. De Craene and G. Berx. Regulatory networks defining emt during cancer initiationand progression.
Nature Reviews Cancer , 13(2):97, 2013.[9] J. Lamouille, S.and Xu and R. Derynck. Molecular mechanisms of epithelial–mesenchymal transition.
Nature Reviews Molecular Cell Biology , 15(3):178, 2014.[10] M. Lu, M. K. Jolly, H. Levine, J. N. Onuchic, and E. Ben-Jacob. MicroRNA-based regu-lation of epithelial–hybrid–mesenchymal fate determination.
Proceedings of the NationalAcademy of Sciences , 110(45):18144–18149, 2013.[11] W. Kolch, M. Halasz, M. Granovskaya, and B. N. Kholodenko. The dynamic control ofsignal transduction networks in cancer cells.
Nature Reviews Cancer , 15(9):515, 2015.[12] M. Hirsch. Differential equations and convergence almost everywhere in strongly mono-tone flows.
Contemporary Mathematics , 17:267–285, 1983.[13] M. Hirsch. Systems of differential equations that are competitive or cooperative ii:Convergence almost everywhere.
SIAM J. Mathematical Analysis , 16:423–439, 1985.[14] M. Hirsch. Stability and convergence in strongly monotone dynamical systems.
J. ReineAngew. Math. , 383:1–53, 1988.[15] M. Hirsch and H.L. Smith.
Monotone dynamical systems . Elsevier, Amsterdam, 2005.Handbook of Differential Equations, Ordinary Differential Equations: Volume 2.[16] H. L. Smith.
Monotone dynamical systems: An introduction to the theory of competi-tive and cooperative systems . AMS, Providence, RI, 1995. Mathematical Surveys andMonographs, vol. 41.[17] D. Angeli, J. E. Ferrell, and E. D. Sontag. Detection of multistability, bifurcations, andhysteresis in a large class of biological positive-feedback systems.
Proceedings of theNational Academy of Sciences , 101(7):1822–1827, 2004.[18] E. D. Sontag. Monotone and near-monotone biochemical networks.
Systems and Syn-thetic biology , 1(2):59–87, 2007. 1419] D. Angeli, P. De Leenheer, and E. D. Sontag. Graph-theoretic characterizations ofmonotonicity of chemical networks in reaction coordinates.
Journal of MathematicalBiology , 61(4):581–616, 2010.[20] L. V. Sharova, A. A. Sharov, T. Nedorezov, Y. Piao, N. Shaik, and Minoru S. H. Ko.Database for mRNA half-life of 19 977 genes obtained by DNA microarray analysis ofpluripotent and differentiating mouse embryonic stem cells.
DNA Research , 16(1):45–58,2008.[21] Y. Guo, J. Liu, S. J. Elfenbein, Y. Ma, M. Zhong, C. Qiu, Y. Ding, and J. Lu. Char-acterization of the mammalian miRNA turnover landscape.
Nucleic Acids Research ,43(4):2326–2341, 2015.[22] L. Wang and E.D. Sontag. Singularly perturbed monotone systems and an applicationto double phosphorylation cycles.
J. Nonlinear Sciences , 2008. DOI 10.1007/s00332-008-9021-2.[23] L. Wang and E. D. Sontag. Singularly perturbed monotone systems and an applicationto double phosphorylation cycles.
Journal of Nonlinear Science , 18(5):527–550, 2008.[24] M. Lu, M. K. Jolly, R. Gomoto, B. Huang, J. Onuchic, and E. Ben-Jacob. Tristabilityin cancer-associated microRNA-TF chimera toggle switch.
The Journal of PhysicalChemistry B , 117(42):13164–13174, 2013.[25] P. ´Erdi and J. T´oth.
Mathematical models of chemical reactions: theory and applicationsof deterministic and stochastic models . Manchester University Press, 1989.[26] M. Ali Al-Radhawi, D. Angeli, and E. D. Sontag. A computational framework fora Lyapunov-enabled analysis of biochemical reaction networks.
PLoS ComputationalBiology , 16(2): e1007681, 2020.[27] M. Ali Al-Radhawi, D. Del Vecchio, and E. D. Sontag. Multi-modality in gene regulatorynetworks with slow promoter kinetics.
PLoS Computational Biology , 15(2):e1006784,2019.[28] C. Ender and G. Meister. Argonaute proteins at a glance.
J Cell Sci , 123(11):1819–1823,2010.[29] A. Berman and R. J. Plemmons.
Nonnegative matrices in the mathematical sciences .SIAM, 1994.[30] N. Fenichel. Geometric singular perturbation theory for ordinary differential equations.
J. of Differential Equations , 31:53–98, 1979.[31] C. K. R. T. Jones. Geometric singular perturbation theory. In
Dynamical Systems(Montecatini. Terme), Lect. Notes in Math. 1609 . Springer-Verlag, Berlin, 1994.[32] K. Nipp. Smooth attractive invariant manifolds of singularly perturbed ODEs. In