Dynamical Properties of a Two-gene Network with Hysteresis
arXiv:1402.0130v1 [cs.SY] 1 Feb 2014
Dynamical Properties of a Two-gene Network with Hysteresis
Qin Shu [email protected]
Ricardo G. Sanfelice [email protected]
February 4, 2014 Technical Report Hybrid Dynamics and Control Laboratory Department of Aerospace and Mechanical Engineering University of Arizona, Tucson Technical Report No. UA/AME/HDC-2014-001. Status: NOT PUBLISHED. Readers of this material have the responsibility to inform all of the authors promptly if they wish to reuse, modify, correct, publish, or distribute any portion of this report. http://www.u.arizona.edu/∼sricardo/index.php?n=Main.TechnicalReports
Dynamical Properties of a Two-gene Network with Hysteresis Qin Shu and Ricardo G. Sanfelice∗ February 4, 2014
Abstract A mathematical model for a two-gene regulatory network is derived and several of their properties analyzed. Due to the presence of mixed continuous/discrete dynamics and hysteresis, we employ a hybrid systems model to capture the dynamics of the system. The proposed model incorporates binary hysteresis with different thresholds capturing the interaction between the genes. We analyze properties of the solutions and asymptotic stability of equilibria in the system as a function of its parameters. Our analysis reveals the presence of limit cycles for a certain range of parameters, behavior that is associated with hysteresis. The set of points defining the limit cycle is characterized and its asymptotic stability properties are studied. Furthermore, the stability property of the limit cycle is robust to small perturbations. Numerical simulations are presented to illustrate the results.
∗
Q. Shu and R. G. Sanfelice are with the Department of Aerospace and Mechanical Engineering, University of Arizona 1130 N. Mountain Ave, AZ 85721. Email: [email protected], [email protected]. This research has been partially supported by the National Science Foundation under CAREER Grant no. ECS-1150306 and by the Air Force Office of Scientific Research under Grant no. FA9550-12-1-0366.
1
Contents 1 Introduction 1.1 Mathematical modeling of genetic regulatory networks . . . . . . . . . . . . 1.2 The role of hysteresis in genetic regulatory networks . . . . . . . . . . . . . . 1.3 Contributions and organization of the paper . . . . . . . . . . . . . . . . . .
3 3 3 3
2 A Hybrid Model for Genetic Networks w/Hysteresis 2.1 Introduction to Hybrid System Modeling . . . . . . . . . . . . . . . . . . . . 2.2 Modeling of a Two-Gene Network . . . . . . . . . . . . . . . . . . . . . . . .
4 6 7
3 Dynamical Properties of the Two-Gene Hybrid System Model 3.1 Existence of solutions . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Characterization of equilibria . . . . . . . . . . . . . . . . . . . . 3.3 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Asymptotic stability of isolated equilibrium points . . . . . 3.3.2 Stability properties of the limit cycle . . . . . . . . . . . . 3.4 Robustness properties . . . . . . . . . . . . . . . . . . . . . . . . 4 Numerical results 4.1 Isolated equilibrium points 4.1.1 Case 1 of Table 1 . 4.1.2 Case 2 of Table 1 . 4.1.3 Case 3 of Table 1 . 4.1.4 Case 4 of Table 1 . 4.2 Equilibrium set S . . . . .
in Table . . . . . . . . . . . . . . . . . . . . . . . . .
1 . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . .
10 10 11 15 15 15 20
. . . . . .
21 21 21 22 23 24 25
5 Conclusion
31
A
34 34 40 44 52
A.1 A.2 A.3 A.4
Proof Proof Proof Proof
of of of of
Proposition Proposition Proposition Proposition
3.3 3.1 3.5 3.6
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
2
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
1 1.1
Introduction Mathematical modeling of genetic regulatory networks
In recent years, the development of advanced experimental techniques in molecular biology has led to a growing interest in mathematical modeling methods for the study of genetic regulatory networks; see [1] for a literature review. A number of gene regulatory network models have been proposed to capture their main properties [2], [3], [4], [5], [6], [7], [8]. Boolean models capture the dynamics of the discrete switch in genetic networks. As introduced by Glass and Kauffman in [3], Boolean regulation functions, typically modeled as sigmoidal or step functions, can be combined with linear system models to enforce certain logic rules. The properties of such a class of piecewise linear models have been studied in the mathematical biology literature, e.g., [4, 5, 2, 6]. Snoussi presented a discrete mapping approach in [4] to study the qualitative properties of the dynamics of genetic regulatory networks. In this work, the properties of the discrete mapping were studied to determine stable isolated steady states as well as limit cycles. In [5], Gouz´ e and Sari employ the concept of Filippov solution to study piecewise linear models of genetic regulatory networks with discontinuities occurring on hyperplanes defined by thresholds on the variables. Chaves and coauthors [2] studied the robustness of Boolean models of gene control networks. de Jong and coauthors [6] presented a method for qualitative simulation of genetic regulatory networks based on the piecewise linear model of [3]. Genetic regulatory networks with continuous dynamics coupled with switching can be written as a hybrid system. In [7] and [8], the authors apply hybrid systems tools to model a variety of cell biology problems. More recently, hybrid models have been used in [9] for the study of molecular interactions. It is important to note that hysteresis behavior, which is typically present in genetic regulatory networks, has not been considered in the models mentioned above.
1.2
The role of hysteresis in genetic regulatory networks
Hysteresis is an important phenomenon in genetic regulatory networks. It is characterized by behavior in which, for instance, once a gene has been inhibited due to the concentration of cellular protein reaching a particularly low value, a higher value of cellular protein concentration is required to express it. In his survey paper on the impact of genetic modeling on tumorigenesis and drug discovery [10], Huang stated that “hysteresis is a feature that a synthetic model has to capture.” Through experiments, Das and coauthors [11] demonstrated the existence of hysteresis in lymphoid cells and the interaction of continuous evolution of some cellular proteins. Hysteresis was also found to be present in mammalian genetic regulatory networks; see, e.g., [12, 13]. More importantly, it has been observed that hysteresis is a key mechanism contributing to oscillatory behavior in computational biological models [14], [15]. On the other hand, it is well known that hysteresis is one of the key factors that makes a system robust to noise and parametric uncertainties [16], [17].
1.3
Contributions and organization of the paper
Our work is motivated by the following facts: 3
1. Piecewise linear models do not incorporate hysteresis, although it plays a key role in the dynamics of genetic regulatory networks. In fact, as we establish in this paper, hysteresis leads to oscillatory, robust behavior in two-gene networks. 2. The discontinuities introduced by the Boolean regulation functions yield a non-smooth dynamical system, for which classical analysis tools cannot be applied to study existence of solutions, stability, robustness, etc. Motivated by these two limitations, we propose a hybrid system model that captures both continuous and discrete dynamics of genetic regulatory networks with hysteresis behavior. We combine the methodology of piecewise linear modeling of genetic regulatory networks with the framework of hybrid dynamical systems in [18], and construct a hybrid system model for a genetic network with two genes. Our model incorporates hysteresis explicitly, which we found leads to limit cycles. We prove existence of solutions and compute the equilibrium points in terms of parameters for the system. We analyze the stability of the isolated equilibrium points and determined conditions under which a limit cycle exists. It is found that hysteresis is the key mechanism leading to hysteresis, as without hysteresis, the limit cycle converges to an isolated equilibrium point (cf. [4]). The stability of the limit cycle is established using a novel approach consisting of measuring distance between solutions of hybrid systems (rather than the distance to the limit cycle as in classical continuous-time systems). Moreover, we show that the asymptotic stability of the limit cycle is robust to small perturbations. The remainder of this paper is organized as follows. In Section 2, a mathematical framework of hybrid dynamical system is introduced and then applied to model a two-gene network. The analysis of existence of solutions, stability, and robustness are presented in Section 3. Section 4 presents simulations validating our results.
2
A Hybrid Systems Model for Genetic Regulatory Networks with Hysteresis
Models of genetic regulatory networks given by piecewise-linear differential equations have been proposed in [8], [19]. Such models take the form 1 x˙ = f (x) − γx,
x ≥ 0,
(1)
where x = [x1 , x2 , . . . , xn ]⊤ and xi represents the concentration of the protein in the i-th cell, f = [f1 , f2 , . . . , fn ]⊤ is a function, γ = [γ1 , γ2, . . . , γn ]⊤ is a vector of constants, and 1 ≤ i ≤ n. For each i, fi is a function representing the rate of synthesis, while γi represents the degradation rate constant P of the protein. The function fi is typically defined as the linear combination fi (x) = ℓ∈L kiℓ biℓ (x) where kiℓ is the nonzero and nonnegative growth rate constants, biℓ is a Boolean regulation function that describes the gene regulation logic, and L = {1, 2 . . . , n} is the set of indices of regulation functions. The modeling strategy for the Boolean regulation functions bil is a key element that captures the behavior of a particular genetic regulatory network. A major feature of a 1
The notation x ≥ 0 is equivalent to xi ≥ 0 for each i.
4
genetic regulatory network is the presence of threshold-like relationships between the system variables, i.e., if a variable xi is above (or below) a certain level, it could cause little or no effect on another variable xj , whereas if xi is below (or above) this certain value, the effect on xj would become more significant (for example, it may increase the value of xj or inhibit the growth of the value of xj ). Boolean regulation functions can be modeled by sigmoidal or step functions, an approach that was first proposed by Glass and Kauffmann [3]. When modeling as a step function, the functions biℓ are given by the combination (linear or nonlinear) of 1 if xi ≥ θ + s (xi , θ) = , s− (xi , θ) = 1 − s+ (xi , θ), (2) 0 if xi < θ where s+ (xi , θ) represents the logic for gene expression when the protein concentration exceeds a threshold θ, while s− (xi , θ) represents the logic for gene inhibition. To illustrate this modeling approach, let us consider the genetic regulatory network shown in Figure 1. Genes a and b encode proteins A and B , respectively. When the concentration of protein A is below certain threshold, it will inhibit gene b. Similarly, protein B inhibits gene a when the concentration of protein B is above certain threshold. In this way, a set of piecewise-linear differential equations representing the behavior in Figure 1 is given by x˙ 2 = k2 s+ (x1 , θ1 ) − γ2 x2 ,
x˙ 1 = k1 s− (x2 , θ2 ) − γ1 x1 ,
(3)
where x1 is representing the concentration of protein A, while x2 is the concentration of protein B . The constants θ1 , θ2 are the thresholds associated with concentrations of protein A and B, respectively. B A
a
b
Figure 1: A genetic regulatory network of two genes (a and b), each encoding for a protein (A and B). Lines ending in arrows represent genetic expression triggers, while lines ending in flatheads refer to genetic inhibition triggers. In this model, gene a is expressed at a rate k1 when x2 is below the threshold θ2 . Similarly, gene b is expressed at a rate k2 when x1 is above the threshold θ1 . Degradations of both proteins are assumed to be proportional to their own concentrations, a mechanism that is captured by −γ1 x1 and −γ2 x2 , respectively. Note that the model in (3) capturing the interaction between gene a and gene b does not incorporate binary hysteresis. Furthermore, due to the discontinuities introduced by the Boolean regulation functions, it is not straightforward to argue that solutions to (3) exist 5
from every initial value of x. In order to overcome such limitations, we propose a hybrid system with hysteresis for this two gene genetic regulatory network, to which hybrid systems tools for analysis of existence of solutions and asymptotic stability can be applied.
2.1
Introduction to Hybrid System Modeling
Following [18] and [20], a hybrid system in this paper is defined by four objects: • A set C ⊂ Rn , called the flow set. • A set D ⊂ Rn , called the jump set. • A single-valued mapping F : Rn → Rn , called the flow map. • A set-valued mapping G: Rn ⇒ Rn , called the jump map. The flow map F defines the continuous dynamics on the flow set C, while the jump map G defines the discrete dynamics or jumps on the jump set D. These objects are referred to as the data of the hybrid system H. Then, defining z ∈ Rn to be the state of the system, H can be written in the compact form z˙ = F (z) z ∈ C H: z + ∈ G(z) z ∈ D Solutions to hybrid systems are given by hybrid arcs which are trajectories defined on hybrid time domains. Definition 2.1 (hybrid time domain) A set E is a hybrid time domain if for all (T, J) ∈ E, E ∩ ([0, T ] × {0, 1, ..., J}) is a compact hybrid time domain; that is, it can be written as j−1 ∪j=0 ([tj , tj+1 ], j) for some finite sequence of times 0 ≤ t0 ≤ t1 ≤ . . . ≤ tj . Definition 2.2 (hybrid arc) A hybrid arc φ is a function that takes values from Rn , is defined on a hybrid time domain dom φ, and is such that t 7→ φ(t, j) is locally absolutely continuous for every j, (t, j) ∈ dom φ. Hybrid time domains impose a specific structure on the domains of solutions to hybrid systems. In simple words, solutions to H are defined on intervals of flow [tj , tj+1] indexed by the jump time j when tj+1 > tj . Hybrid arcs specify the functions that define solutions to hybrid systems when the following conditions are satisfied. We refer the reader to [20, 18] for more details on the definition of solutions to hybrid systems. Definition 2.3 (solution) A hybrid arc φ is a solution to the hybrid system H if φ(0, 0) ∈ C ∪ D and (S1) For all j ∈ N := {0, 1, 2, . . .} and almost all t such that (t, j) ∈ dom φ, φ(t, j) ∈ C,
φ(t, j) = F (φ(t, j))
(S2) For all (t, j) ∈ dom φ such that (t, j + 1) ∈ dom φ, φ(t, j) ∈ D,
φ(t, j + 1) ∈ G(φ(t, j)) 6
Solutions to hybrid systems are classified as follows: • A solution φ to H is said to be nontrivial if dom φ contains at least two points. • A solution φ to H is said to be complete if dom φ is unbounded. • A solution φ to H is said to be Zeno if it is complete and the projection of dom φ onto Rn≥0 is bounded. • A solution φ to H is said to be maximal if there does not exist another solution ϕ to H such that dom ϕ is a proper subset of dom φ, and ϕ(t, j) = φ(t, j) for all (t, j) ∈ dom φ. The reader is referred to [18] and [20] for more details on this hybrid system framework.
2.2
Modeling of a Two-Gene Network
To model the genetic network in (3) as a hybrid system H, two discrete logic variables, q1 and q2 , are introduced. The dynamics of these variables depend on the thresholds, θ1 and θ2 , respectively. As one of our goals is to introduce binary hysteresis in the model in (3), we define hysteresis level constants h1 and h2 associated with gene a and gene b, respectively. In this way, qi is governed by dynamics such that the evolution in Figure 2 holds. qi 1
0
xi θi − hi
θi + hi
θimax
Figure 2: The update mechanism of qi as a function of xi and previous values of qi . The state of the hybrid system is defined as z = [x1 , x2 , q1 , q2 ]⊤ , where z ∈ Z := R2≥0 ×{0, 1}2 ; x1 , x2 are (nonnegative) continuous states representing protein concentrations; and q1 , q2 are discrete variables. Here, R≥0 := [0, +∞). We specify constants θ1 and θ2 , usually inferred from biological data, satisfying 0 < θ1 < θ1max , 0 < θ2 < θ2max , where θ1max and θ2max are the maximal value of the concentration of protein A and of the protein B , respectively. To define the continuous dynamics of the hybrid system capturing the evolution of (3), we rewrite the piecewise-linear differential equation (3) by replacing the s+ term with the 7
logic variables qi , and the s− term with the complement of the logic variable qi , i.e., 1 − qi . Note that the discrete logic variables qi only change at jumps, i.e., they are constants during flows. Then, q˙i = 0. In this way, the continuous dynamics are governed by the differential equation x˙ 1 = k1 (1 − q2 ) − γ1 x1 , x˙ 2 = k2 q1 − γ2 x2 , q˙1 = q˙2 = 0, from where we obtain the flow map k1 (1 − q2 ) − γ1 x1 k2 q1 − γ2 x2 . F (z) = 0 0
(4)
Now, we describe the discrete update of the state vector z, i.e., we define G and D. To illustrate this construction, we explain how to model the mechanism in Figure 2 for q1 . When q1 = 0 and x1 = θ1 + h1 the state q1 is updated to 1. We write this update law as q1+ = 1. When q1 = 1
x1 = θ1 − h1 ,
and
then the state q1 is updated to 0, i.e., q1+ = 0. It follows that the mechanism of q1 in Figure 2 can be captured by triggering jumps when the components of z satisfy q1 = 0, x1 = θ1 + h1
or
q1 = 1, x1 = θ1 − h1
Note that the update mechanism for q2 is similar to that of q1 just discussed. We can define the flow and jump sets in a compact form by defining functions η1 (x1 , q1 ) := (2q1 − 1)(−x1 + θ1 + (1 − 2q1 )h1 ) η2 (x2 , q2 ) := (2q2 − 1)(−x2 + θ2 + (1 − 2q2 )h2 ). In this way, the flow set is given by C := {z ∈ Z : η1 (x1 , q1 ) ≤ 0, η2 (x2 , q2 ) ≤ 0}
(5)
and the jump set is given by D = {z ∈ C : η1 (x1 , q1 ) = 0} ∪ {z ∈ C : η2 (x2 , q2 ) = 0} 8
(6)
To define the jump map, first note that at jumps, the continuous states x1 and x2 do not change. Then, we conveniently define x1 x1 x2 , g2 (z) := x2 , g1 (z) := q1 1 − q1 1 − q2 q2
so that the jump map G is given by g1 (z) G(z) := g2 (z) {g1 (z), g2 (z)}
if η1 (x1 , q1 ) = 0, η2 (x2 , q2 ) < 0 if η1 (x1 , q1 ) < 0, η2 (x2 , q2 ) = 0 if η1 (x1 , q1 ) = 0, η2 (x2 , q2 ) = 0.
The above definitions determine a hybrid system for (3), which is given by k (1 − q ) − γ x 1 2 1 1 k2 q1 − γ2 x2 z∈C z˙ = F (z) = 0 H:z∈Z 0 + z ∈ G(z) z ∈ D,
(7)
(8)
where C is in (5), G is in (7), and D is in (6). Its parameters are given by the positive constants k1 , k2 , γ1 , γ2 , θ1 , θ2 , h1 , h2 , which satisfy θ1 + h1 < θ1max , θ2 + h2 < θ2max , θ1 − h1 > 0, θ2 − h2 > 0. Figure 3 depicts a hybrid automaton representation of this system when sequentially transitioning between (q1 , q2 ) = (0, 0), (1, 0), (1, 1), (0, 1). Lemma 2.4 The data (C, F, D, G) satisfies the following conditions: (A1) The sets C and D are closed. (A2) The map z 7→ F (z) is continuous on C. (A3) The set-valued mapping z 7→ G(z) is outer semicontinuous2 relative to R4 and locally bounded, and, for all z ∈ D, G(z) is nonempty. Proof: Properties (A1) and (A2) are obvious. Property (A3) holds since the graph of G, which is given by {(x, y) : y ∈ G(z) } , is closed. 2
A set-valued mapping G : S ⇒ Rn with S ⊂ Rn is outer semicontinuous relative to S if for any z ∈ S ∞ and any sequence {zi }∞ i=1 with zi ∈ S, limi→∞ zi = z, and any sequence {wi }i=1 with wi ∈ G(zi ) and limi→∞ wi = w we have w ∈ G(z).
9
x1 = θ1 − h1 , x2 > θ2 − h2 x˙ 1 = −γ1 x1 x˙ 2 = −γ2 x2 q1 = 0 q2 = 1
x˙ 1 = −γ1 x1 x˙ 2 = k2 − γ2 x2 q1 = 1 q2 = 1
x1 < θ1 + h1 x2 = θ2 − h2
x1 > θ1 − h1 x2 = θ2 + h2
x˙ 1 = k1 − γ1 x1 x˙ 2 = −γ2 x2 q1 = 0 q2 = 0
x˙ 1 = k1 − γ1 x1 x˙ 2 = k2 − γ2 x2 q1 = 1 q2 = 0
x1 = θ1 + h1 , x2 < θ2 + h2
Figure 3: A hybrid automaton representation of the two-gene genetic regulatory network for sequential transitions of (q1 , q2 ).
3
Dynamical Properties of the Two-Gene Hybrid System Model
3.1
Existence of solutions
Proposition 3.1 From every point in C ∪D, there exists a nontrivial solution for the hybrid system H in (8). Furthermore, every maximal solution is complete and the projection of its hybrid time domain on R≥0 is unbounded, i.e., every solution is not Zeno. The proof of this result uses the conditions for the existence of solutions to H in [18] for general hybrid systems. More precisely, consider the hybrid system H and let z(0, 0) ∈ C ∪ D. If z(0, 0) ∈ D or (VC) there exists a neighborhood U of z(0, 0) such that3 for every z ∈ U ∩ C, F (z) ∩ TC (z) 6= ∅, then there exists a nontrivial solution to H from z(0, 0). If (VC) holds for every z(0, 0) ∈ C \ D, then there exists a nontrivial solution to H from every initial point in C ∪ D, and every maximal solution z satisfies exactly one of the following conditions: 3
TC (z) denotes the tangent cone of C at z, i.e., it is the set of all v for which there exists a sequence of real numbers αi ց 0 and a sequence vi → v such that for every i = 1, 2, ..., x + αi vi ∈ C.
10
1. z is complete; 2. dom z is bounded and the last interval is of the form [tJ , tJ+1 ), where J = sup(t,j)∈dom z j has nonempty interior and t 7→ φ(t, J) is a maximal solution to z˙ = F (z), in fact limt→T |z(t, J)|=∞, where T = sup(t,j)∈dom z t; 3. z(T, J) ∈ / C ∪ D, where (T, J) = sup dom z. Furthermore, if G(D) ⊂ C ∪ D, then 3) above does not occur. The proof of Proposition 3.1 uses these conditions and is given in A.2.
3.2
Characterization of equilibria
We compute the set of isolated equilibrium points z ∗ as well as (nonisolated, dense) sets of equilibria for the hybrid system H in (8). For general hybrid systems, isolated equilibrium points are points that are an isolated equilibrium point of z˙ ∈ F (z), z ∈ C or of z + ∈ G(z), z ∈ D. On the other hand, an equilibrium set (not necessarily an isolated equilibrium point) for a hybrid system H is defined as a set that is (strongly) forward invariant. Definition 3.2 (Equilibrium set) A set S ⊂ C ∪ D is an equilibrium set of H if for every initial condition z(0, 0) ∈ S, every solution z to H satisfies z(t, j) ∈ S for all (t, j) ∈ S. The following results determine the equilibria of (8) for a range of parameters of the system. Proposition 3.3 The equilibria of the hybrid system H in (8) is given in Table 1 in terms of the positive constants k1 , k2 , γ1 , γ2 , θ1 , θ1max , θ2 , θ2max , h1 , and h2 satisfying the conditions therein. The set S ⊂ C ∪ D in case 5 is an equilibrium set and is given by S=
4 [
Si ,
i=1
where4 S1 :=
S2
4
x ∈ R2 :
:= x ∈ R2 : (
S3 := S4
(
x ∈ R2 :
:= x ∈ R2 :
# ) − p0 (1) exp(−γ1 s) , s ∈ [0, t′1 ] × {(0, 0)} x= p0 (2) exp(−γ2 s) k1 k1 − − p (1) exp(−γ s) 1 1 γ1 γ1 ′ x= , s ∈ [0, t2 ] × {(1, 0)} k2 k2 − − p (2) exp(−γ s) 1 2 γ2 γ2 # ) " p2 (1) exp(−γ 1 s) ′ , s ∈ [0, t3 ] × {(1, 1)} x = k2 − kγ22 − p2 (2) exp(−γ2 s) γ2 p3 (1) exp(−γ1 s) ′ x= , s ∈ [0, t4 ] × {(0, 1)} p3 (2) exp(−γ2 s) "
k1 γ1
−
k1 γ1
pi (j) is the j-th component of pi .
11
(9)
and p0 , p1 , p2 , p3 ∈ R2 are the vertices of the set S(see Figure 4), where t′1
"
k1 γ1
− p0 (1)
# γ1
"
1
, t′2
k2 γ2
− p1 (2)
= ln k2 − (θ1 + h1 ) − (θ2 + h2 ) γ2 1 1 p2 (1) γ1 ′ p3 (2) γ2 ′ t3 = ln , t4 = ln , θ1 − h1 θ2 − h2 = ln
k1 γ1
# γ1
2
,
and p0 =
"
γγ1 # 2 2 (θ1 − h1 ) θp23−h (2) , θ2 − h2
(10)
θ1 + h1 γ2 k1 p1 = −(θ1 +h1 ) γ1 , γ1 (θ2 − h2 ) k1
γ1
p2 = p3 =
k1 γ1
−
"
k2 γ2
k1 γ1
(11)
−p0 (1)
k2 −(θ2 +h2 ) γγ12 γ2 − (θ1 + h1 ) k2 , −p1 (2) γ2
(12)
θ2 + h2
−
θ1 − h1 k2 γ2
− (θ2 + h2 )
Moreover, the period of the limit cycle is given by
θ1 −h1 p2 (1)
γ2 γ1
#
.
T = t′1 + t′2 + t′3 + t′4 .
(13)
(14)
The following result provides a more constructive characterization of S. Corollary 3.4 Under the conditions of Proposition 3.3, if furthermore, γ1 = γ2 = γ, then p0 (1) =
−d6 + d8 + d7 − d5 − d4 − d3 + d2 , d1
where d1 = 2h1 k22 γ + h2 k1 k2 γ + k1 k2 γθ2 − 2h1 h2 k2 γ 2 −2h1 k2 γ 2 θ2 d2 = k1 k2 γθ1 θ2 , d3 = h2 k1 k2 γθ1 , d4 = h1 k1 k2 γθ2 , d5 = h1 h2 k1 k2 γ, d7 = h2 k12 k2 , d8 = h1 k1 k22
12
(15)
Table 1: Equilibria of the hybrid system (8).
1 2 3 4 5
Conditions on constants θ1 + h1 < kγ11 < θ1max 0 < γk22 < θ2 + h2 0
0,γ2 > 0, then f4 (x) points outside of the C4 .
−γ1 x1 . Since −γ2 (θ2 − h2 )
Then, since there is no isolated equilibrium point in C4 , the trajectories starting in C4 have x component that flow towards {0 ≤ x1 ≤ θ1 + h1 , x2 = θ2 − h2 }. 37
Combining the above arguments, Figure 18 shows the transition sequence of q ∈ Q for case 5 in Table 1. Now, we compute the value of the trajectories as they transition according to the said sequence. The differential equation for the x components of the continuous dynamics of H can be evaluated for each possible value of q and written as x˙ = K(q) − Γx, where Γ =
γ1 0 0 γ2
(34)
and K : Q → R2×1 is given by k 1 0 0 0 K(q) = k1 k2 0 k2
if q = (0, 0), if q = (0, 1), if q = (1, 0), if q = (1, 1).
Restricted to C, (34) is a linear time-invariant system. For any initial condition z(0, 0) = [x(0, 0)⊤ q(0, 0)⊤ ]⊤ ∈ C, the unique solution to (34) for each t ≥ 0, up to the first jump, is given by x(t, 0) = u(q(0, 0)) + exp(−Γt)(x(0, 0) − u(q(0, 0))), (35) where u(q) = Γ−1 K(q) for each q ∈ Q. 01
11
00
10
Figure 18: Transition graph of q ∈ Q. 0 and the initial value of x given by For the initial value of q given by q 0 = 0 p0 (1) , x(0, 0) = p0 = θ2 − h2 the solution to (34) is given by x(t, 0) =
u1 − (u1 − x1 (0, 0)) exp(−γ1 t) x2 (0, 0) exp(−γ2 t) 38
k1 k1 γ1 Since K((0, 0)) = implies that u(q) = . Note that x2 (t, 0) where u1 = 0 0 decreases to zero in C1 and that x1 (t, 0) is increasing in C1 , reaching the threshold value h i1 1 −x1 (0,0) γ1 . A jump of q to occurs at t = t1 , j = 0. x1 = θ1 + h1 at t′1 = ln uu11−(θ 1 +h1 ) 0 After the jump, the initial value of x is p1 , where " # θ1 + h1 γ γ2 . p1 = (36) 1 1 +h1 ) p0 (2) u1u−(θ 1 −p0 (1) k1 . γ1
Proceeding similarly as when the initial state was p0 , we obtain the following expressions for p2 , p3 , and p4 : " γ1 # u2 −p2 (2) γ2 u − (u − p (1)) 1 1 1 u2 −p1 (2) p2 = , θ2 + h2 # " " (37) γ1 # θ1 − h1 p4 (2) γ2 γ p (1) 2 , p4 = 3 p3 (2) , p3 = (1) γ1 u2 − (u2 − p2 (2)) pp23 (1) θ2 − h2
where u2 =
k2 . γ2
Also, similarly, we have the expression for t′2 , t′3 , t′4 :
t′2 = ln
h
u2 −p1 (2) u2 −(θ2 +h2 )
i γ1
2
,
t′3 = ln
h
p2 (1) θ1 −h1
i γ1
1
,
t′4 = ln
h
p3 (2) θ2 −h2
i γ1
2
.
Then, the period of the limit cycle is given by T = t′1 + t′2 + t′3 + t′4 . Note that t1 = t′1 , t2 = t′1 + t′2 , t3 = t′1 + t′2 + t′3 and t4 = t′1 + t′2 + t′3 + t′4 , where t1 , t2 , t3 and t4 define the jump times (t1 , 0), (t2 , 1), (t3 , 2), and (t4 , 3). Now, we define the map ρ : [0, θ1max ] → R as ρ(r) = ρ4 ◦ ρ3 ◦ ρ2 ◦ ρ1 (r), where ρ1 (r) = (θ2 − h2 )
u1 −(θ1 +h1 ) u1 −r
γγ2
ρ2 (r) = u1 − (u1 − (θ1 + h1 ))
Then, r such that
ρ3 (r) = u2 − (u2 − (θ2 + h2 )) γ1 2 γ2 ρ4 (r) = (θ1 − h1 ) θ2 −h . r
1
,
u2 −(θ2 +h2 ) u2 −r
γ2 θ1 −h1 γ 1
r
,
γγ1
2
,
ρ(r) = r
defines p0 (1). Then, combining the above expressions, we obtain (10)-(13). Finally, using (35), the set S is constructed by combining the x components of the (unique) solutions between these points. Since each piece of the x component corresponding to a constant value of q is a solution to a linear system, this set of points has the property that, from every point in it, the only existing solution from that point stays in the set, i.e., the set is strongly forward invariant. 39
A.2
Proof of Proposition 3.1
To verify the sufficient conditions for the existence of nontrivial solutions from an initial point in C ∪ D, it is enough to show that F (z) ∈ TC (z) for every z ∈ C \ D in the boundary of C (the (VC) condition holds for every point in the interior of C.) Next, we consider each possible case. 1. Let z ∈ C1 \ D1 = {z : q1 = 0, q2 = 0, 0 ≤ x1 < θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }. Let TC11 (z) TC21 (z) TC31 (z) TC41 (z) TC51 (z) TC61 (z) TC71 (z) TC81 (z)
= = = = = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w1 : w1 : w1 : w2 : w1 : w1 : w1 : w2
≥ 0, w2 ≥ 0}, ≥ 0, w2 ≥ 0}, ≤ 0, w2 ≤ 0}, ≤ 0, w2 ≤ 0}.
≤ 0}, ≥ 0}, ≥ 0}, ≤ 0},
Then, the tangent cone of C1 at points z = (x, q) is given as follows: • For x ∈ {x : x1 = 0, x2 = θ2 + h2 }, TC1 (z) = TC11 (z), • For x ∈ {x : x1 = 0, 0 < x2 < θ2 + h2 }, TC1 (z) = TC21 (z). • For x ∈ {x : x1 = 0, x2 = 0}, TC1 (z) = TC31 (z). • For x ∈ {0 < x1 < θ1 + h1 , x2 = 0}, TC1 (z) = TC41 (z), • For x ∈ {x : x1 = θ1 + h1 , x2 = 0}, TC1 (z) = TC51 (z). • For x ∈ {x1 = θ1 + h1 , 0 < x2 < θ2 + h2 }, TC1 (z) = TC61 (z), • For x ∈ {x1 = θ1 + h1 , x2 = θ2 + h2 }, TC1 (z) = TC71 (z), • For x ∈ {0 < x1 < θ1 + h1 , x2 = θ2 + h2 }, TC1 (z) = TC81 (z). Now, we check the vector field F on the boundary of C1 away from D1 . k1 • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) := . Since k1 > 0, −γ2 x2 f1 (x) points inside of C1 .
40
x2
q1 = 0
q2 = 0
θ2max C4 8
1
θ2 + h2
C3 7
θ2 − h2 2
6
C1
3
4 θ1 − h1
C2
5 θ1 + h1
x1 θ1max
Figure 19: Tangent cones on the boundaries of C1 . • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = 0}, f1 (x) := tangent to the boundary of C1 .
k1 − γ1 x1 . Then, f1 (x) is 0
Then, F (z) ∈ TC1 holds, implying that (VC) holds at each point z ∈ C1 \ D1 . 2. Let z ∈ C2 \ D2 = {z : q1 = 1, q2 = 0, x1 ≥ θ1 − h1 , 0 ≤ x2 < θ2 + h2 }. Let TC12 (z) TC22 (z) TC32 (z) TC42 (z) TC52 (z)
= = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w2 : w1 : w1 : w1 : w2
≥ 0}, ≥ 0, w2 ≥ 0}, ≥ 0}, ≤ 0, w2 ≥ 0}, ≤ 0}.
Then, the tangent cone of C2 is given by as follows: • For x ∈ {x : x1 > θ1 − h1 , x2 = 0}, TC2 (z) = TC12 (z), • For x ∈ {x : x1 = θ1 − h1 , x2 = 0}, TC2 (z) = TC22 (z). • For x ∈ {x : x1 = θ1 − h1 , 0 ≤ x2 ≤ θ2 + h2 }, TC2 (z) = TC32 (z). • For x ∈ {x1 = θ1 − h1 , x2 = θ2 + h2 }, TC2 (z) = TC42 (z), • For x ∈ {x : x1 > θ1 − h1 , x2 = θ2 + h2 }, TC2 (z) = TC52 (z).
41
Figure 20 depicts the tangent cones on the boundaries of C2 . x2 q1 = 1 q2 = 0 θ2max
θ2 + h2
1
5
2
C2
θ2 − h2
3 θ1−h1
θ1+h1
x1
4
θ1max
Figure 20: The tangent cone of C2 when θ1 + h1
θ1 − h1 , x2 > θ2 − h2 }. Since there are no points in the boundary of C3 that are not in D3 , (VC) holds for free. 4. Let z ∈ C4 \ D4 = {z : q1 = 0, q2 = 1, 0 ≤ x1 ≤ θ1 + h1 , x2 > θ2 − h2 }. Let TC14 (z) TC24 (z) TC34 (z) TC44 (z) TC54 (z)
= = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w1 : w1 : w2 : w1 : w1
≤ 0}, ≤ 0, w2 ≥ 0}, ≥ 0}, ≥ 0, w2 ≥ 0}, ≥ 0}.
Then, the tangent cone of C4 is given by as follows: • For x ∈ {x : x1 = θ1 + h1 , x2 > θ2 − h2 }, TC4 (z) = TC14 (z), • For x ∈ {x : x1 = θ1 + h1 , x2 = θ2 − h2 }, TC4 (z) = TC24 (z).
42
• For x ∈ {x : 0 < x1 < θ1 + h1 , x2 = θ2 − h2 }, TC4 (z) = TC34 (z). • For x ∈ {x1 = 0, x2 = θ2 − h2 }, TC4 (z) = TC44 (z), • For x ∈ {x : x1 = 0, x2 > θ2 − h2 }, TC4 (z) = TC54 (z). x2 q2= 1
q1= 0 θ2max
5
C4
θ2 −h2 4
3
1
θ2 +h2 2
x1 θ1 −h1 θ1 +h1
θ1max
Figure 21: Tangent cones on the boundaries of C4 . Now, we check the vector field F on the boundary of C4 away from D4 . 0 , and f4 (x) is tangent • When x ∈ {x : x1 = 0, x2 ≥ θ2 − h2 }, f4 (x) := −γ2 x2 to the boundary of C4 . Then, F (z) ∈ TC4 (z) holds at each point z ∈ C4 \ D4 . Combining the above arguments, for each case in Table 1, (VC) holds and nontrivial solutions to H in (8) exist. Since fi is linear, every solution to z˙ = F (z) subject to z ∈ C does not escape to infinity by flowing. Then, condition 2) below (VC) does not hold. Since G(D) ⊂ C ∪ D, then condition 3) therein does not hold either. Finally, every solution is not Zeno, since at most after the second jump, every solution needs to flow (linearly) from the value after the jump, which is given by G, for at least 2 min{h1 , h2 } in the x1 or in the x2 direction (for certain initial conditions, e.g., z(0, 0) = [θ1 − h1 , θ2 − h2 , 1, 1], solutions jump twice consecutively, and after that, flow for the said amount).
43
A.3
Proof of Proposition 3.5
Note that when θ1 + h1
θ2 + h2 . γ2
0 = x˙ 1 = k1 − γ1 x1 0 = x˙ 2 = γ2 x2 ⇒
k1 γ1
θ1max
< θ1 + h1 . (a) Case
⇒ k1 − γ1 x1 = 0 ⇒ x∗1 = −γ2 x2 = 0 ⇒ x∗2 = 0,
k1 γ1
from where we have z2∗ = [x∗ ⊤ 0 0]⊤ ∈ C1 . Now, change to e coordinates given by e1 = x1 − x∗1 ,
e2 = x2 − x∗2 .
We have that e˙ 1 = = e˙ 2 = =
x˙ 1 − x˙ ∗1 = k1 − γ1 x1 − 0 = k1 − γ1 (e1 + x∗1 ) k1 − γ1 e1 − k1 = −γ1 e1 , x˙ 2 = k2 − γ2 x2 = k2 − γ2 (e2 + x∗2 ) k2 − γ2 e2 − k2 = −γ2 e2 .
Then, the x component of z˙ = F (z) becomes −γ1 0 e˙ = e. 0 −γ2 . Then, since γ1 and γ2 are positive, we have that z2∗ is stable (this property can be certified with the Lyapunov function V (e) = e⊤ e). Now we check the vector fields on the boundaries of C1 (see Figure 26). k1 − γ1 (θ1 + h1 ) . Since • When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 0 < γk11 < θ1 − h1 , we have that f1 (x) points inside C1 .
47
x2 q1 = 1
q2 = 1
θ2max C4
C3
C1
C2
θ2 + h2 θ2 − h2
x1 θ1−h1
θ1+h1
θ1max
Figure 26: The vector field at the boundaries of C1 , when 0
0, we have • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 that f1 (x) points inside C1 . 0
θ1 − h1 , x2 = θ2 + h2 }, f2 (x) =
– if
48
k1 − γ1 x1 . Since 0 < kγ11 < θ1 − h1 , we have that f2 (x) points inside of k2 − γ2 (θ2 + h2 ) C2 . Since there is no equilibrium point in C2 for this range of parameters and the dynamics of x are linear, x1 reaches θ1 − h1 for every point z ∈ C2 . Then, a jump occurs. After the jump, the solution belongs to C1 .
k2 γ2
> θ2 + h2 (see Figure 25(b)), when x ∈ {x : x1 = θ1 − h1 , 0 ≤ x2 ≤ k1 − γ1 (θ1 − h1 ) θ2 + h2 }, f2 (x) = . Since 0 < kγ11 < θ1 − h1 , we have that k2 − γ2 x2 f2 (x) points outside of C2 . When x ∈ {x : x1 ≥ θ1 − h1 , x2 = θ2 + h2 }, f1 (x) = k1 − γ1 x1 . Since 0 < γk11 < θ1 − h1 , we have that f2 (x) points outside k2 − γ2 (θ2 + h2 ) C2 .
– If
Then, from z(0, 0) ∈ C2 , solutions will leave C2 and jump into C1 or C3 . • For every initial point z(0, 0) ∈ C3 : k2 γ2
< θ2 −h2 , when x ∈ {x : x1 ≥ θ1 −h1 , x2 = θ2 −h2 } (see similar case shown in −γ1 x1 . We have that f3 (x) points outside of Figure 24(a)), f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , C3 . When x ∈ {x : x1 = θ1 −h1 , x2 ≥ θ2 −h2 }, we have f3 (x) = k2 − γ2 x2 which points outside of C3 .
– if
k2 γ2
< θ2 + h2 (see similar case shown in Figure 24(b)), when x ∈ {x : −γ1 (θ1 − h1 ) points outside of C3 . When x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, f3 (x) = k2 − γ2 x2 −γ1 x1 . We have that x ∈ {x : x1 > θ1 − h1 , x2 = θ2 − h2 }, f3 (x) = k2 − γ2 (θ2 − h2 ) f3 (x) points inside of C3 .
– If θ2 − h2
θ2 − h2 }, f4 (x) = −γ2 x2 f4 (x) points inside C4 . Then, for every initial condition z(0, 0) ∈ C4 , solutions will reach C1 in finite time. From the above analysis, we have that: 1) from C2 trajectories go to either C1 or C3 ; 2) from C3 trajectories go to either C1 or C4 ; 3) from C4 trajectories go to C1 . Then, trajectories eventually enter C1 , which, using the arguments above, implies that the equilibrium point z2∗ is globally asymptotically stable. 49
Similarly, for case 4 in Table 1, it can be proven that z2∗ is globally asymptotically stable. Since stability of z2∗ was proven in the proof of case 2 in Table 1, now we check the vector fields on the boundaries of C1 when θ1 − h1 < γk11 < θ1 + h1 , θ2 + h2 < γk22 < θ2max (see Figure 27). x2 q1 = 0
q2 = 0
θ2max
C4
C3
C1
C2
θ2 + h2 θ2 − h2
x1 θ1−h1 θ1+h1
θ1max
Figure 27: The vector field at the boundaries of C1 , when θ1 − h1 < k2 < θ2max . γ2
• When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = θ1 − h1
0, we have that f2 (x) points inside of C2 . Then, once a trajectory enters or starts from C2 , it will stay or never leave C2 Then, the set C2 is forward invariant. Since the equilibrium point z1∗ belongs to C2 , every trajectory reaching or starting from C2 converges to z1∗ . For z(0, 0) ∈ C1 , we check the vector fields on the boundaries of C1 (see Figure 31(b)). k1 − γ1 (θ1 + h1 ) . Since • When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 θ1 − h1 < γk11 < θ1 + h1 , we have that f1 (x) points inside C1 . k1 − γ1 x1 . Since • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = θ2 + h2 }, f1 (x) = −γ2 (θ2 + h2 ) θ1 − h1 < γk11 < θ1 + h1 , we have that f1 (x) points inside C1 . k1 − γ1 (θ1 + h1 ) . f1 (x) is • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = 0}, f1 (x) = 0 tangent to the boundary of C1 . k1 . Since k1 > 0, we have • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 that f1 (x) points inside C1 . 53
x2
x2
q1 = 1
q1 = 0
q2 = 0
q2 = 0
θ2max
θ2max C4
C3
C4
C3
C1
C2
θ2 + h2 θ2 − h2
θ2 + h2 θ2 − h2 C2
C1
x1
x1 θ1−h1 θ1+h1 (a)
θ1−h1 θ1+h1 (b)
θ1max
θ1max
Figure 31: when θ1 −h1 < kγ11 < θ1 +h1 ,0 < γk22 < θ2 +h2 : (a) the vector field at the boundaries of C2 , (b) the vector field at the boundaries of C1 . Then, once a trajectory enters or starts from C1 , it will stay or never leave C1 . Then, the set C1 is forward invariant. Since the equilibrium point z2∗ belongs to C1 , every trajectory reaching or starting from C1 converges to z2∗ . For initial points in C4 , the vector field of the boundary at C4 is shownin Figure 23 (b). −γ1 x1 . Then, we have When x ∈ {x : 0 < x1 < θ1 + h1 , x2 = θ2 − h2 }, f4 (x) = −γ2 (θ2 − h2 ) that f4 (x) points outside C4 and every solution leaves C4 by jumping into C1 . Then, for every initial condition z(0, 0) ∈ C4 , solutions will reach C1 in finite time. As z2∗ ∈ C1 , so the trajectory will stay in C1 and converge to z2∗ . If z(0, 0) ∈ C3 , the vector field at the boundary of C3 is similar as that shown in Figure 24. 1. If the parameters are in the range of θ1 − h1 < kγ11 < θ1 + h1 , 0 < γk22 < θ2 − h2 (see max similar case shown in Figure 24(a)), when x ∈ {x : θ1 − h1 < x1 < θ1 , x2 = θ2 − h2 }, −γ1 x1 . We have that f3 (x) points outside of C3 . When x ∈ {x : f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , which points outside x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, we have f3 (x) = k2 − γ2 x2 C3 . Depending on which jump set the solution hits, two possibilities of the equilibrium points exist. • If the trajectory hits the set x ∈ {x : x1 = θ1 − h1 , x2 > θ2 − h2 }, which leads to an update of q1 , the solution will jump into C4 . For this case, the trajectory will reach C1 in finite time and converge to z2∗ . • If the trajectory hits the set x ∈ {x : x1 > θ1 − h1 , x2 = θ2 − h2 }, q2 is updated to 0 from 1, the trajectory will enter C2 , and converge to the equilibrium point z1∗ , which is inside of C2 . 54
2. If the parameters are in the range of θ1 − h1 < kγ11 < θ1 + h1 , θ2 − h2 < kγ22 < θ2 + h2 (see similar 24(b)), when x ∈ {x : x1 > θ1 + h1 , x2 = θ2 − h2 }, case shown in Figure −γ1 x1 . We have that f3 (x) points inside of C3 . When x ∈ {x : f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , which points outside x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, we have f3 (x) = k2 − γ2 x2 of C3 . Thus, when the trajectory hits the set {x : x1 = θ1 − h1 , x2 > θ2 − h2 }, a jump will occur, the solution will jump into C4 and finally enter C1 and converge to z2∗ .
55
Dynamical Properties of a Two-gene Network with Hysteresis
Qin Shu [email protected]
Ricardo G. Sanfelice [email protected]
February 4, 2014 Technical Report Hybrid Dynamics and Control Laboratory Department of Aerospace and Mechanical Engineering University of Arizona, Tucson Technical Report No. UA/AME/HDC-2014-001. Status: NOT PUBLISHED. Readers of this material have the responsibility to inform all of the authors promptly if they wish to reuse, modify, correct, publish, or distribute any portion of this report. http://www.u.arizona.edu/∼sricardo/index.php?n=Main.TechnicalReports
Dynamical Properties of a Two-gene Network with Hysteresis Qin Shu and Ricardo G. Sanfelice∗ February 4, 2014
Abstract A mathematical model for a two-gene regulatory network is derived and several of their properties analyzed. Due to the presence of mixed continuous/discrete dynamics and hysteresis, we employ a hybrid systems model to capture the dynamics of the system. The proposed model incorporates binary hysteresis with different thresholds capturing the interaction between the genes. We analyze properties of the solutions and asymptotic stability of equilibria in the system as a function of its parameters. Our analysis reveals the presence of limit cycles for a certain range of parameters, behavior that is associated with hysteresis. The set of points defining the limit cycle is characterized and its asymptotic stability properties are studied. Furthermore, the stability property of the limit cycle is robust to small perturbations. Numerical simulations are presented to illustrate the results.
∗
Q. Shu and R. G. Sanfelice are with the Department of Aerospace and Mechanical Engineering, University of Arizona 1130 N. Mountain Ave, AZ 85721. Email: [email protected], [email protected]. This research has been partially supported by the National Science Foundation under CAREER Grant no. ECS-1150306 and by the Air Force Office of Scientific Research under Grant no. FA9550-12-1-0366.
1
Contents 1 Introduction 1.1 Mathematical modeling of genetic regulatory networks . . . . . . . . . . . . 1.2 The role of hysteresis in genetic regulatory networks . . . . . . . . . . . . . . 1.3 Contributions and organization of the paper . . . . . . . . . . . . . . . . . .
3 3 3 3
2 A Hybrid Model for Genetic Networks w/Hysteresis 2.1 Introduction to Hybrid System Modeling . . . . . . . . . . . . . . . . . . . . 2.2 Modeling of a Two-Gene Network . . . . . . . . . . . . . . . . . . . . . . . .
4 6 7
3 Dynamical Properties of the Two-Gene Hybrid System Model 3.1 Existence of solutions . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Characterization of equilibria . . . . . . . . . . . . . . . . . . . . 3.3 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Asymptotic stability of isolated equilibrium points . . . . . 3.3.2 Stability properties of the limit cycle . . . . . . . . . . . . 3.4 Robustness properties . . . . . . . . . . . . . . . . . . . . . . . . 4 Numerical results 4.1 Isolated equilibrium points 4.1.1 Case 1 of Table 1 . 4.1.2 Case 2 of Table 1 . 4.1.3 Case 3 of Table 1 . 4.1.4 Case 4 of Table 1 . 4.2 Equilibrium set S . . . . .
in Table . . . . . . . . . . . . . . . . . . . . . . . . .
1 . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . .
10 10 11 15 15 15 20
. . . . . .
21 21 21 22 23 24 25
5 Conclusion
31
A
34 34 40 44 52
A.1 A.2 A.3 A.4
Proof Proof Proof Proof
of of of of
Proposition Proposition Proposition Proposition
3.3 3.1 3.5 3.6
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
2
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
1 1.1
Introduction Mathematical modeling of genetic regulatory networks
In recent years, the development of advanced experimental techniques in molecular biology has led to a growing interest in mathematical modeling methods for the study of genetic regulatory networks; see [1] for a literature review. A number of gene regulatory network models have been proposed to capture their main properties [2], [3], [4], [5], [6], [7], [8]. Boolean models capture the dynamics of the discrete switch in genetic networks. As introduced by Glass and Kauffman in [3], Boolean regulation functions, typically modeled as sigmoidal or step functions, can be combined with linear system models to enforce certain logic rules. The properties of such a class of piecewise linear models have been studied in the mathematical biology literature, e.g., [4, 5, 2, 6]. Snoussi presented a discrete mapping approach in [4] to study the qualitative properties of the dynamics of genetic regulatory networks. In this work, the properties of the discrete mapping were studied to determine stable isolated steady states as well as limit cycles. In [5], Gouz´ e and Sari employ the concept of Filippov solution to study piecewise linear models of genetic regulatory networks with discontinuities occurring on hyperplanes defined by thresholds on the variables. Chaves and coauthors [2] studied the robustness of Boolean models of gene control networks. de Jong and coauthors [6] presented a method for qualitative simulation of genetic regulatory networks based on the piecewise linear model of [3]. Genetic regulatory networks with continuous dynamics coupled with switching can be written as a hybrid system. In [7] and [8], the authors apply hybrid systems tools to model a variety of cell biology problems. More recently, hybrid models have been used in [9] for the study of molecular interactions. It is important to note that hysteresis behavior, which is typically present in genetic regulatory networks, has not been considered in the models mentioned above.
1.2
The role of hysteresis in genetic regulatory networks
Hysteresis is an important phenomenon in genetic regulatory networks. It is characterized by behavior in which, for instance, once a gene has been inhibited due to the concentration of cellular protein reaching a particularly low value, a higher value of cellular protein concentration is required to express it. In his survey paper on the impact of genetic modeling on tumorigenesis and drug discovery [10], Huang stated that “hysteresis is a feature that a synthetic model has to capture.” Through experiments, Das and coauthors [11] demonstrated the existence of hysteresis in lymphoid cells and the interaction of continuous evolution of some cellular proteins. Hysteresis was also found to be present in mammalian genetic regulatory networks; see, e.g., [12, 13]. More importantly, it has been observed that hysteresis is a key mechanism contributing to oscillatory behavior in computational biological models [14], [15]. On the other hand, it is well known that hysteresis is one of the key factors that makes a system robust to noise and parametric uncertainties [16], [17].
1.3
Contributions and organization of the paper
Our work is motivated by the following facts: 3
1. Piecewise linear models do not incorporate hysteresis, although it plays a key role in the dynamics of genetic regulatory networks. In fact, as we establish in this paper, hysteresis leads to oscillatory, robust behavior in two-gene networks. 2. The discontinuities introduced by the Boolean regulation functions yield a non-smooth dynamical system, for which classical analysis tools cannot be applied to study existence of solutions, stability, robustness, etc. Motivated by these two limitations, we propose a hybrid system model that captures both continuous and discrete dynamics of genetic regulatory networks with hysteresis behavior. We combine the methodology of piecewise linear modeling of genetic regulatory networks with the framework of hybrid dynamical systems in [18], and construct a hybrid system model for a genetic network with two genes. Our model incorporates hysteresis explicitly, which we found leads to limit cycles. We prove existence of solutions and compute the equilibrium points in terms of parameters for the system. We analyze the stability of the isolated equilibrium points and determined conditions under which a limit cycle exists. It is found that hysteresis is the key mechanism leading to hysteresis, as without hysteresis, the limit cycle converges to an isolated equilibrium point (cf. [4]). The stability of the limit cycle is established using a novel approach consisting of measuring distance between solutions of hybrid systems (rather than the distance to the limit cycle as in classical continuous-time systems). Moreover, we show that the asymptotic stability of the limit cycle is robust to small perturbations. The remainder of this paper is organized as follows. In Section 2, a mathematical framework of hybrid dynamical system is introduced and then applied to model a two-gene network. The analysis of existence of solutions, stability, and robustness are presented in Section 3. Section 4 presents simulations validating our results.
2
A Hybrid Systems Model for Genetic Regulatory Networks with Hysteresis
Models of genetic regulatory networks given by piecewise-linear differential equations have been proposed in [8], [19]. Such models take the form 1 x˙ = f (x) − γx,
x ≥ 0,
(1)
where x = [x1 , x2 , . . . , xn ]⊤ and xi represents the concentration of the protein in the i-th cell, f = [f1 , f2 , . . . , fn ]⊤ is a function, γ = [γ1 , γ2, . . . , γn ]⊤ is a vector of constants, and 1 ≤ i ≤ n. For each i, fi is a function representing the rate of synthesis, while γi represents the degradation rate constant P of the protein. The function fi is typically defined as the linear combination fi (x) = ℓ∈L kiℓ biℓ (x) where kiℓ is the nonzero and nonnegative growth rate constants, biℓ is a Boolean regulation function that describes the gene regulation logic, and L = {1, 2 . . . , n} is the set of indices of regulation functions. The modeling strategy for the Boolean regulation functions bil is a key element that captures the behavior of a particular genetic regulatory network. A major feature of a 1
The notation x ≥ 0 is equivalent to xi ≥ 0 for each i.
4
genetic regulatory network is the presence of threshold-like relationships between the system variables, i.e., if a variable xi is above (or below) a certain level, it could cause little or no effect on another variable xj , whereas if xi is below (or above) this certain value, the effect on xj would become more significant (for example, it may increase the value of xj or inhibit the growth of the value of xj ). Boolean regulation functions can be modeled by sigmoidal or step functions, an approach that was first proposed by Glass and Kauffmann [3]. When modeling as a step function, the functions biℓ are given by the combination (linear or nonlinear) of 1 if xi ≥ θ + s (xi , θ) = , s− (xi , θ) = 1 − s+ (xi , θ), (2) 0 if xi < θ where s+ (xi , θ) represents the logic for gene expression when the protein concentration exceeds a threshold θ, while s− (xi , θ) represents the logic for gene inhibition. To illustrate this modeling approach, let us consider the genetic regulatory network shown in Figure 1. Genes a and b encode proteins A and B , respectively. When the concentration of protein A is below certain threshold, it will inhibit gene b. Similarly, protein B inhibits gene a when the concentration of protein B is above certain threshold. In this way, a set of piecewise-linear differential equations representing the behavior in Figure 1 is given by x˙ 2 = k2 s+ (x1 , θ1 ) − γ2 x2 ,
x˙ 1 = k1 s− (x2 , θ2 ) − γ1 x1 ,
(3)
where x1 is representing the concentration of protein A, while x2 is the concentration of protein B . The constants θ1 , θ2 are the thresholds associated with concentrations of protein A and B, respectively. B A
a
b
Figure 1: A genetic regulatory network of two genes (a and b), each encoding for a protein (A and B). Lines ending in arrows represent genetic expression triggers, while lines ending in flatheads refer to genetic inhibition triggers. In this model, gene a is expressed at a rate k1 when x2 is below the threshold θ2 . Similarly, gene b is expressed at a rate k2 when x1 is above the threshold θ1 . Degradations of both proteins are assumed to be proportional to their own concentrations, a mechanism that is captured by −γ1 x1 and −γ2 x2 , respectively. Note that the model in (3) capturing the interaction between gene a and gene b does not incorporate binary hysteresis. Furthermore, due to the discontinuities introduced by the Boolean regulation functions, it is not straightforward to argue that solutions to (3) exist 5
from every initial value of x. In order to overcome such limitations, we propose a hybrid system with hysteresis for this two gene genetic regulatory network, to which hybrid systems tools for analysis of existence of solutions and asymptotic stability can be applied.
2.1
Introduction to Hybrid System Modeling
Following [18] and [20], a hybrid system in this paper is defined by four objects: • A set C ⊂ Rn , called the flow set. • A set D ⊂ Rn , called the jump set. • A single-valued mapping F : Rn → Rn , called the flow map. • A set-valued mapping G: Rn ⇒ Rn , called the jump map. The flow map F defines the continuous dynamics on the flow set C, while the jump map G defines the discrete dynamics or jumps on the jump set D. These objects are referred to as the data of the hybrid system H. Then, defining z ∈ Rn to be the state of the system, H can be written in the compact form z˙ = F (z) z ∈ C H: z + ∈ G(z) z ∈ D Solutions to hybrid systems are given by hybrid arcs which are trajectories defined on hybrid time domains. Definition 2.1 (hybrid time domain) A set E is a hybrid time domain if for all (T, J) ∈ E, E ∩ ([0, T ] × {0, 1, ..., J}) is a compact hybrid time domain; that is, it can be written as j−1 ∪j=0 ([tj , tj+1 ], j) for some finite sequence of times 0 ≤ t0 ≤ t1 ≤ . . . ≤ tj . Definition 2.2 (hybrid arc) A hybrid arc φ is a function that takes values from Rn , is defined on a hybrid time domain dom φ, and is such that t 7→ φ(t, j) is locally absolutely continuous for every j, (t, j) ∈ dom φ. Hybrid time domains impose a specific structure on the domains of solutions to hybrid systems. In simple words, solutions to H are defined on intervals of flow [tj , tj+1] indexed by the jump time j when tj+1 > tj . Hybrid arcs specify the functions that define solutions to hybrid systems when the following conditions are satisfied. We refer the reader to [20, 18] for more details on the definition of solutions to hybrid systems. Definition 2.3 (solution) A hybrid arc φ is a solution to the hybrid system H if φ(0, 0) ∈ C ∪ D and (S1) For all j ∈ N := {0, 1, 2, . . .} and almost all t such that (t, j) ∈ dom φ, φ(t, j) ∈ C,
φ(t, j) = F (φ(t, j))
(S2) For all (t, j) ∈ dom φ such that (t, j + 1) ∈ dom φ, φ(t, j) ∈ D,
φ(t, j + 1) ∈ G(φ(t, j)) 6
Solutions to hybrid systems are classified as follows: • A solution φ to H is said to be nontrivial if dom φ contains at least two points. • A solution φ to H is said to be complete if dom φ is unbounded. • A solution φ to H is said to be Zeno if it is complete and the projection of dom φ onto Rn≥0 is bounded. • A solution φ to H is said to be maximal if there does not exist another solution ϕ to H such that dom ϕ is a proper subset of dom φ, and ϕ(t, j) = φ(t, j) for all (t, j) ∈ dom φ. The reader is referred to [18] and [20] for more details on this hybrid system framework.
2.2
Modeling of a Two-Gene Network
To model the genetic network in (3) as a hybrid system H, two discrete logic variables, q1 and q2 , are introduced. The dynamics of these variables depend on the thresholds, θ1 and θ2 , respectively. As one of our goals is to introduce binary hysteresis in the model in (3), we define hysteresis level constants h1 and h2 associated with gene a and gene b, respectively. In this way, qi is governed by dynamics such that the evolution in Figure 2 holds. qi 1
0
xi θi − hi
θi + hi
θimax
Figure 2: The update mechanism of qi as a function of xi and previous values of qi . The state of the hybrid system is defined as z = [x1 , x2 , q1 , q2 ]⊤ , where z ∈ Z := R2≥0 ×{0, 1}2 ; x1 , x2 are (nonnegative) continuous states representing protein concentrations; and q1 , q2 are discrete variables. Here, R≥0 := [0, +∞). We specify constants θ1 and θ2 , usually inferred from biological data, satisfying 0 < θ1 < θ1max , 0 < θ2 < θ2max , where θ1max and θ2max are the maximal value of the concentration of protein A and of the protein B , respectively. To define the continuous dynamics of the hybrid system capturing the evolution of (3), we rewrite the piecewise-linear differential equation (3) by replacing the s+ term with the 7
logic variables qi , and the s− term with the complement of the logic variable qi , i.e., 1 − qi . Note that the discrete logic variables qi only change at jumps, i.e., they are constants during flows. Then, q˙i = 0. In this way, the continuous dynamics are governed by the differential equation x˙ 1 = k1 (1 − q2 ) − γ1 x1 , x˙ 2 = k2 q1 − γ2 x2 , q˙1 = q˙2 = 0, from where we obtain the flow map k1 (1 − q2 ) − γ1 x1 k2 q1 − γ2 x2 . F (z) = 0 0
(4)
Now, we describe the discrete update of the state vector z, i.e., we define G and D. To illustrate this construction, we explain how to model the mechanism in Figure 2 for q1 . When q1 = 0 and x1 = θ1 + h1 the state q1 is updated to 1. We write this update law as q1+ = 1. When q1 = 1
x1 = θ1 − h1 ,
and
then the state q1 is updated to 0, i.e., q1+ = 0. It follows that the mechanism of q1 in Figure 2 can be captured by triggering jumps when the components of z satisfy q1 = 0, x1 = θ1 + h1
or
q1 = 1, x1 = θ1 − h1
Note that the update mechanism for q2 is similar to that of q1 just discussed. We can define the flow and jump sets in a compact form by defining functions η1 (x1 , q1 ) := (2q1 − 1)(−x1 + θ1 + (1 − 2q1 )h1 ) η2 (x2 , q2 ) := (2q2 − 1)(−x2 + θ2 + (1 − 2q2 )h2 ). In this way, the flow set is given by C := {z ∈ Z : η1 (x1 , q1 ) ≤ 0, η2 (x2 , q2 ) ≤ 0}
(5)
and the jump set is given by D = {z ∈ C : η1 (x1 , q1 ) = 0} ∪ {z ∈ C : η2 (x2 , q2 ) = 0} 8
(6)
To define the jump map, first note that at jumps, the continuous states x1 and x2 do not change. Then, we conveniently define x1 x1 x2 , g2 (z) := x2 , g1 (z) := q1 1 − q1 1 − q2 q2
so that the jump map G is given by g1 (z) G(z) := g2 (z) {g1 (z), g2 (z)}
if η1 (x1 , q1 ) = 0, η2 (x2 , q2 ) < 0 if η1 (x1 , q1 ) < 0, η2 (x2 , q2 ) = 0 if η1 (x1 , q1 ) = 0, η2 (x2 , q2 ) = 0.
The above definitions determine a hybrid system for (3), which is given by k (1 − q ) − γ x 1 2 1 1 k2 q1 − γ2 x2 z∈C z˙ = F (z) = 0 H:z∈Z 0 + z ∈ G(z) z ∈ D,
(7)
(8)
where C is in (5), G is in (7), and D is in (6). Its parameters are given by the positive constants k1 , k2 , γ1 , γ2 , θ1 , θ2 , h1 , h2 , which satisfy θ1 + h1 < θ1max , θ2 + h2 < θ2max , θ1 − h1 > 0, θ2 − h2 > 0. Figure 3 depicts a hybrid automaton representation of this system when sequentially transitioning between (q1 , q2 ) = (0, 0), (1, 0), (1, 1), (0, 1). Lemma 2.4 The data (C, F, D, G) satisfies the following conditions: (A1) The sets C and D are closed. (A2) The map z 7→ F (z) is continuous on C. (A3) The set-valued mapping z 7→ G(z) is outer semicontinuous2 relative to R4 and locally bounded, and, for all z ∈ D, G(z) is nonempty. Proof: Properties (A1) and (A2) are obvious. Property (A3) holds since the graph of G, which is given by {(x, y) : y ∈ G(z) } , is closed. 2
A set-valued mapping G : S ⇒ Rn with S ⊂ Rn is outer semicontinuous relative to S if for any z ∈ S ∞ and any sequence {zi }∞ i=1 with zi ∈ S, limi→∞ zi = z, and any sequence {wi }i=1 with wi ∈ G(zi ) and limi→∞ wi = w we have w ∈ G(z).
9
x1 = θ1 − h1 , x2 > θ2 − h2 x˙ 1 = −γ1 x1 x˙ 2 = −γ2 x2 q1 = 0 q2 = 1
x˙ 1 = −γ1 x1 x˙ 2 = k2 − γ2 x2 q1 = 1 q2 = 1
x1 < θ1 + h1 x2 = θ2 − h2
x1 > θ1 − h1 x2 = θ2 + h2
x˙ 1 = k1 − γ1 x1 x˙ 2 = −γ2 x2 q1 = 0 q2 = 0
x˙ 1 = k1 − γ1 x1 x˙ 2 = k2 − γ2 x2 q1 = 1 q2 = 0
x1 = θ1 + h1 , x2 < θ2 + h2
Figure 3: A hybrid automaton representation of the two-gene genetic regulatory network for sequential transitions of (q1 , q2 ).
3
Dynamical Properties of the Two-Gene Hybrid System Model
3.1
Existence of solutions
Proposition 3.1 From every point in C ∪D, there exists a nontrivial solution for the hybrid system H in (8). Furthermore, every maximal solution is complete and the projection of its hybrid time domain on R≥0 is unbounded, i.e., every solution is not Zeno. The proof of this result uses the conditions for the existence of solutions to H in [18] for general hybrid systems. More precisely, consider the hybrid system H and let z(0, 0) ∈ C ∪ D. If z(0, 0) ∈ D or (VC) there exists a neighborhood U of z(0, 0) such that3 for every z ∈ U ∩ C, F (z) ∩ TC (z) 6= ∅, then there exists a nontrivial solution to H from z(0, 0). If (VC) holds for every z(0, 0) ∈ C \ D, then there exists a nontrivial solution to H from every initial point in C ∪ D, and every maximal solution z satisfies exactly one of the following conditions: 3
TC (z) denotes the tangent cone of C at z, i.e., it is the set of all v for which there exists a sequence of real numbers αi ց 0 and a sequence vi → v such that for every i = 1, 2, ..., x + αi vi ∈ C.
10
1. z is complete; 2. dom z is bounded and the last interval is of the form [tJ , tJ+1 ), where J = sup(t,j)∈dom z j has nonempty interior and t 7→ φ(t, J) is a maximal solution to z˙ = F (z), in fact limt→T |z(t, J)|=∞, where T = sup(t,j)∈dom z t; 3. z(T, J) ∈ / C ∪ D, where (T, J) = sup dom z. Furthermore, if G(D) ⊂ C ∪ D, then 3) above does not occur. The proof of Proposition 3.1 uses these conditions and is given in A.2.
3.2
Characterization of equilibria
We compute the set of isolated equilibrium points z ∗ as well as (nonisolated, dense) sets of equilibria for the hybrid system H in (8). For general hybrid systems, isolated equilibrium points are points that are an isolated equilibrium point of z˙ ∈ F (z), z ∈ C or of z + ∈ G(z), z ∈ D. On the other hand, an equilibrium set (not necessarily an isolated equilibrium point) for a hybrid system H is defined as a set that is (strongly) forward invariant. Definition 3.2 (Equilibrium set) A set S ⊂ C ∪ D is an equilibrium set of H if for every initial condition z(0, 0) ∈ S, every solution z to H satisfies z(t, j) ∈ S for all (t, j) ∈ S. The following results determine the equilibria of (8) for a range of parameters of the system. Proposition 3.3 The equilibria of the hybrid system H in (8) is given in Table 1 in terms of the positive constants k1 , k2 , γ1 , γ2 , θ1 , θ1max , θ2 , θ2max , h1 , and h2 satisfying the conditions therein. The set S ⊂ C ∪ D in case 5 is an equilibrium set and is given by S=
4 [
Si ,
i=1
where4 S1 :=
S2
4
x ∈ R2 :
:= x ∈ R2 : (
S3 := S4
(
x ∈ R2 :
:= x ∈ R2 :
# ) − p0 (1) exp(−γ1 s) , s ∈ [0, t′1 ] × {(0, 0)} x= p0 (2) exp(−γ2 s) k1 k1 − − p (1) exp(−γ s) 1 1 γ1 γ1 ′ x= , s ∈ [0, t2 ] × {(1, 0)} k2 k2 − − p (2) exp(−γ s) 1 2 γ2 γ2 # ) " p2 (1) exp(−γ 1 s) ′ , s ∈ [0, t3 ] × {(1, 1)} x = k2 − kγ22 − p2 (2) exp(−γ2 s) γ2 p3 (1) exp(−γ1 s) ′ x= , s ∈ [0, t4 ] × {(0, 1)} p3 (2) exp(−γ2 s) "
k1 γ1
−
k1 γ1
pi (j) is the j-th component of pi .
11
(9)
and p0 , p1 , p2 , p3 ∈ R2 are the vertices of the set S(see Figure 4), where t′1
"
k1 γ1
− p0 (1)
# γ1
"
1
, t′2
k2 γ2
− p1 (2)
= ln k2 − (θ1 + h1 ) − (θ2 + h2 ) γ2 1 1 p2 (1) γ1 ′ p3 (2) γ2 ′ t3 = ln , t4 = ln , θ1 − h1 θ2 − h2 = ln
k1 γ1
# γ1
2
,
and p0 =
"
γγ1 # 2 2 (θ1 − h1 ) θp23−h (2) , θ2 − h2
(10)
θ1 + h1 γ2 k1 p1 = −(θ1 +h1 ) γ1 , γ1 (θ2 − h2 ) k1
γ1
p2 = p3 =
k1 γ1
−
"
k2 γ2
k1 γ1
(11)
−p0 (1)
k2 −(θ2 +h2 ) γγ12 γ2 − (θ1 + h1 ) k2 , −p1 (2) γ2
(12)
θ2 + h2
−
θ1 − h1 k2 γ2
− (θ2 + h2 )
Moreover, the period of the limit cycle is given by
θ1 −h1 p2 (1)
γ2 γ1
#
.
T = t′1 + t′2 + t′3 + t′4 .
(13)
(14)
The following result provides a more constructive characterization of S. Corollary 3.4 Under the conditions of Proposition 3.3, if furthermore, γ1 = γ2 = γ, then p0 (1) =
−d6 + d8 + d7 − d5 − d4 − d3 + d2 , d1
where d1 = 2h1 k22 γ + h2 k1 k2 γ + k1 k2 γθ2 − 2h1 h2 k2 γ 2 −2h1 k2 γ 2 θ2 d2 = k1 k2 γθ1 θ2 , d3 = h2 k1 k2 γθ1 , d4 = h1 k1 k2 γθ2 , d5 = h1 h2 k1 k2 γ, d7 = h2 k12 k2 , d8 = h1 k1 k22
12
(15)
Table 1: Equilibria of the hybrid system (8).
1 2 3 4 5
Conditions on constants θ1 + h1 < kγ11 < θ1max 0 < γk22 < θ2 + h2 0
0,γ2 > 0, then f4 (x) points outside of the C4 .
−γ1 x1 . Since −γ2 (θ2 − h2 )
Then, since there is no isolated equilibrium point in C4 , the trajectories starting in C4 have x component that flow towards {0 ≤ x1 ≤ θ1 + h1 , x2 = θ2 − h2 }. 37
Combining the above arguments, Figure 18 shows the transition sequence of q ∈ Q for case 5 in Table 1. Now, we compute the value of the trajectories as they transition according to the said sequence. The differential equation for the x components of the continuous dynamics of H can be evaluated for each possible value of q and written as x˙ = K(q) − Γx, where Γ =
γ1 0 0 γ2
(34)
and K : Q → R2×1 is given by k 1 0 0 0 K(q) = k1 k2 0 k2
if q = (0, 0), if q = (0, 1), if q = (1, 0), if q = (1, 1).
Restricted to C, (34) is a linear time-invariant system. For any initial condition z(0, 0) = [x(0, 0)⊤ q(0, 0)⊤ ]⊤ ∈ C, the unique solution to (34) for each t ≥ 0, up to the first jump, is given by x(t, 0) = u(q(0, 0)) + exp(−Γt)(x(0, 0) − u(q(0, 0))), (35) where u(q) = Γ−1 K(q) for each q ∈ Q. 01
11
00
10
Figure 18: Transition graph of q ∈ Q. 0 and the initial value of x given by For the initial value of q given by q 0 = 0 p0 (1) , x(0, 0) = p0 = θ2 − h2 the solution to (34) is given by x(t, 0) =
u1 − (u1 − x1 (0, 0)) exp(−γ1 t) x2 (0, 0) exp(−γ2 t) 38
k1 k1 γ1 Since K((0, 0)) = implies that u(q) = . Note that x2 (t, 0) where u1 = 0 0 decreases to zero in C1 and that x1 (t, 0) is increasing in C1 , reaching the threshold value h i1 1 −x1 (0,0) γ1 . A jump of q to occurs at t = t1 , j = 0. x1 = θ1 + h1 at t′1 = ln uu11−(θ 1 +h1 ) 0 After the jump, the initial value of x is p1 , where " # θ1 + h1 γ γ2 . p1 = (36) 1 1 +h1 ) p0 (2) u1u−(θ 1 −p0 (1) k1 . γ1
Proceeding similarly as when the initial state was p0 , we obtain the following expressions for p2 , p3 , and p4 : " γ1 # u2 −p2 (2) γ2 u − (u − p (1)) 1 1 1 u2 −p1 (2) p2 = , θ2 + h2 # " " (37) γ1 # θ1 − h1 p4 (2) γ2 γ p (1) 2 , p4 = 3 p3 (2) , p3 = (1) γ1 u2 − (u2 − p2 (2)) pp23 (1) θ2 − h2
where u2 =
k2 . γ2
Also, similarly, we have the expression for t′2 , t′3 , t′4 :
t′2 = ln
h
u2 −p1 (2) u2 −(θ2 +h2 )
i γ1
2
,
t′3 = ln
h
p2 (1) θ1 −h1
i γ1
1
,
t′4 = ln
h
p3 (2) θ2 −h2
i γ1
2
.
Then, the period of the limit cycle is given by T = t′1 + t′2 + t′3 + t′4 . Note that t1 = t′1 , t2 = t′1 + t′2 , t3 = t′1 + t′2 + t′3 and t4 = t′1 + t′2 + t′3 + t′4 , where t1 , t2 , t3 and t4 define the jump times (t1 , 0), (t2 , 1), (t3 , 2), and (t4 , 3). Now, we define the map ρ : [0, θ1max ] → R as ρ(r) = ρ4 ◦ ρ3 ◦ ρ2 ◦ ρ1 (r), where ρ1 (r) = (θ2 − h2 )
u1 −(θ1 +h1 ) u1 −r
γγ2
ρ2 (r) = u1 − (u1 − (θ1 + h1 ))
Then, r such that
ρ3 (r) = u2 − (u2 − (θ2 + h2 )) γ1 2 γ2 ρ4 (r) = (θ1 − h1 ) θ2 −h . r
1
,
u2 −(θ2 +h2 ) u2 −r
γ2 θ1 −h1 γ 1
r
,
γγ1
2
,
ρ(r) = r
defines p0 (1). Then, combining the above expressions, we obtain (10)-(13). Finally, using (35), the set S is constructed by combining the x components of the (unique) solutions between these points. Since each piece of the x component corresponding to a constant value of q is a solution to a linear system, this set of points has the property that, from every point in it, the only existing solution from that point stays in the set, i.e., the set is strongly forward invariant. 39
A.2
Proof of Proposition 3.1
To verify the sufficient conditions for the existence of nontrivial solutions from an initial point in C ∪ D, it is enough to show that F (z) ∈ TC (z) for every z ∈ C \ D in the boundary of C (the (VC) condition holds for every point in the interior of C.) Next, we consider each possible case. 1. Let z ∈ C1 \ D1 = {z : q1 = 0, q2 = 0, 0 ≤ x1 < θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }. Let TC11 (z) TC21 (z) TC31 (z) TC41 (z) TC51 (z) TC61 (z) TC71 (z) TC81 (z)
= = = = = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w1 : w1 : w1 : w2 : w1 : w1 : w1 : w2
≥ 0, w2 ≥ 0}, ≥ 0, w2 ≥ 0}, ≤ 0, w2 ≤ 0}, ≤ 0, w2 ≤ 0}.
≤ 0}, ≥ 0}, ≥ 0}, ≤ 0},
Then, the tangent cone of C1 at points z = (x, q) is given as follows: • For x ∈ {x : x1 = 0, x2 = θ2 + h2 }, TC1 (z) = TC11 (z), • For x ∈ {x : x1 = 0, 0 < x2 < θ2 + h2 }, TC1 (z) = TC21 (z). • For x ∈ {x : x1 = 0, x2 = 0}, TC1 (z) = TC31 (z). • For x ∈ {0 < x1 < θ1 + h1 , x2 = 0}, TC1 (z) = TC41 (z), • For x ∈ {x : x1 = θ1 + h1 , x2 = 0}, TC1 (z) = TC51 (z). • For x ∈ {x1 = θ1 + h1 , 0 < x2 < θ2 + h2 }, TC1 (z) = TC61 (z), • For x ∈ {x1 = θ1 + h1 , x2 = θ2 + h2 }, TC1 (z) = TC71 (z), • For x ∈ {0 < x1 < θ1 + h1 , x2 = θ2 + h2 }, TC1 (z) = TC81 (z). Now, we check the vector field F on the boundary of C1 away from D1 . k1 • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) := . Since k1 > 0, −γ2 x2 f1 (x) points inside of C1 .
40
x2
q1 = 0
q2 = 0
θ2max C4 8
1
θ2 + h2
C3 7
θ2 − h2 2
6
C1
3
4 θ1 − h1
C2
5 θ1 + h1
x1 θ1max
Figure 19: Tangent cones on the boundaries of C1 . • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = 0}, f1 (x) := tangent to the boundary of C1 .
k1 − γ1 x1 . Then, f1 (x) is 0
Then, F (z) ∈ TC1 holds, implying that (VC) holds at each point z ∈ C1 \ D1 . 2. Let z ∈ C2 \ D2 = {z : q1 = 1, q2 = 0, x1 ≥ θ1 − h1 , 0 ≤ x2 < θ2 + h2 }. Let TC12 (z) TC22 (z) TC32 (z) TC42 (z) TC52 (z)
= = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w2 : w1 : w1 : w1 : w2
≥ 0}, ≥ 0, w2 ≥ 0}, ≥ 0}, ≤ 0, w2 ≥ 0}, ≤ 0}.
Then, the tangent cone of C2 is given by as follows: • For x ∈ {x : x1 > θ1 − h1 , x2 = 0}, TC2 (z) = TC12 (z), • For x ∈ {x : x1 = θ1 − h1 , x2 = 0}, TC2 (z) = TC22 (z). • For x ∈ {x : x1 = θ1 − h1 , 0 ≤ x2 ≤ θ2 + h2 }, TC2 (z) = TC32 (z). • For x ∈ {x1 = θ1 − h1 , x2 = θ2 + h2 }, TC2 (z) = TC42 (z), • For x ∈ {x : x1 > θ1 − h1 , x2 = θ2 + h2 }, TC2 (z) = TC52 (z).
41
Figure 20 depicts the tangent cones on the boundaries of C2 . x2 q1 = 1 q2 = 0 θ2max
θ2 + h2
1
5
2
C2
θ2 − h2
3 θ1−h1
θ1+h1
x1
4
θ1max
Figure 20: The tangent cone of C2 when θ1 + h1
θ1 − h1 , x2 > θ2 − h2 }. Since there are no points in the boundary of C3 that are not in D3 , (VC) holds for free. 4. Let z ∈ C4 \ D4 = {z : q1 = 0, q2 = 1, 0 ≤ x1 ≤ θ1 + h1 , x2 > θ2 − h2 }. Let TC14 (z) TC24 (z) TC34 (z) TC44 (z) TC54 (z)
= = = = =
{w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2 {w ∈ R2
: w1 : w1 : w2 : w1 : w1
≤ 0}, ≤ 0, w2 ≥ 0}, ≥ 0}, ≥ 0, w2 ≥ 0}, ≥ 0}.
Then, the tangent cone of C4 is given by as follows: • For x ∈ {x : x1 = θ1 + h1 , x2 > θ2 − h2 }, TC4 (z) = TC14 (z), • For x ∈ {x : x1 = θ1 + h1 , x2 = θ2 − h2 }, TC4 (z) = TC24 (z).
42
• For x ∈ {x : 0 < x1 < θ1 + h1 , x2 = θ2 − h2 }, TC4 (z) = TC34 (z). • For x ∈ {x1 = 0, x2 = θ2 − h2 }, TC4 (z) = TC44 (z), • For x ∈ {x : x1 = 0, x2 > θ2 − h2 }, TC4 (z) = TC54 (z). x2 q2= 1
q1= 0 θ2max
5
C4
θ2 −h2 4
3
1
θ2 +h2 2
x1 θ1 −h1 θ1 +h1
θ1max
Figure 21: Tangent cones on the boundaries of C4 . Now, we check the vector field F on the boundary of C4 away from D4 . 0 , and f4 (x) is tangent • When x ∈ {x : x1 = 0, x2 ≥ θ2 − h2 }, f4 (x) := −γ2 x2 to the boundary of C4 . Then, F (z) ∈ TC4 (z) holds at each point z ∈ C4 \ D4 . Combining the above arguments, for each case in Table 1, (VC) holds and nontrivial solutions to H in (8) exist. Since fi is linear, every solution to z˙ = F (z) subject to z ∈ C does not escape to infinity by flowing. Then, condition 2) below (VC) does not hold. Since G(D) ⊂ C ∪ D, then condition 3) therein does not hold either. Finally, every solution is not Zeno, since at most after the second jump, every solution needs to flow (linearly) from the value after the jump, which is given by G, for at least 2 min{h1 , h2 } in the x1 or in the x2 direction (for certain initial conditions, e.g., z(0, 0) = [θ1 − h1 , θ2 − h2 , 1, 1], solutions jump twice consecutively, and after that, flow for the said amount).
43
A.3
Proof of Proposition 3.5
Note that when θ1 + h1
θ2 + h2 . γ2
0 = x˙ 1 = k1 − γ1 x1 0 = x˙ 2 = γ2 x2 ⇒
k1 γ1
θ1max
< θ1 + h1 . (a) Case
⇒ k1 − γ1 x1 = 0 ⇒ x∗1 = −γ2 x2 = 0 ⇒ x∗2 = 0,
k1 γ1
from where we have z2∗ = [x∗ ⊤ 0 0]⊤ ∈ C1 . Now, change to e coordinates given by e1 = x1 − x∗1 ,
e2 = x2 − x∗2 .
We have that e˙ 1 = = e˙ 2 = =
x˙ 1 − x˙ ∗1 = k1 − γ1 x1 − 0 = k1 − γ1 (e1 + x∗1 ) k1 − γ1 e1 − k1 = −γ1 e1 , x˙ 2 = k2 − γ2 x2 = k2 − γ2 (e2 + x∗2 ) k2 − γ2 e2 − k2 = −γ2 e2 .
Then, the x component of z˙ = F (z) becomes −γ1 0 e˙ = e. 0 −γ2 . Then, since γ1 and γ2 are positive, we have that z2∗ is stable (this property can be certified with the Lyapunov function V (e) = e⊤ e). Now we check the vector fields on the boundaries of C1 (see Figure 26). k1 − γ1 (θ1 + h1 ) . Since • When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 0 < γk11 < θ1 − h1 , we have that f1 (x) points inside C1 .
47
x2 q1 = 1
q2 = 1
θ2max C4
C3
C1
C2
θ2 + h2 θ2 − h2
x1 θ1−h1
θ1+h1
θ1max
Figure 26: The vector field at the boundaries of C1 , when 0
0, we have • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 that f1 (x) points inside C1 . 0
θ1 − h1 , x2 = θ2 + h2 }, f2 (x) =
– if
48
k1 − γ1 x1 . Since 0 < kγ11 < θ1 − h1 , we have that f2 (x) points inside of k2 − γ2 (θ2 + h2 ) C2 . Since there is no equilibrium point in C2 for this range of parameters and the dynamics of x are linear, x1 reaches θ1 − h1 for every point z ∈ C2 . Then, a jump occurs. After the jump, the solution belongs to C1 .
k2 γ2
> θ2 + h2 (see Figure 25(b)), when x ∈ {x : x1 = θ1 − h1 , 0 ≤ x2 ≤ k1 − γ1 (θ1 − h1 ) θ2 + h2 }, f2 (x) = . Since 0 < kγ11 < θ1 − h1 , we have that k2 − γ2 x2 f2 (x) points outside of C2 . When x ∈ {x : x1 ≥ θ1 − h1 , x2 = θ2 + h2 }, f1 (x) = k1 − γ1 x1 . Since 0 < γk11 < θ1 − h1 , we have that f2 (x) points outside k2 − γ2 (θ2 + h2 ) C2 .
– If
Then, from z(0, 0) ∈ C2 , solutions will leave C2 and jump into C1 or C3 . • For every initial point z(0, 0) ∈ C3 : k2 γ2
< θ2 −h2 , when x ∈ {x : x1 ≥ θ1 −h1 , x2 = θ2 −h2 } (see similar case shown in −γ1 x1 . We have that f3 (x) points outside of Figure 24(a)), f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , C3 . When x ∈ {x : x1 = θ1 −h1 , x2 ≥ θ2 −h2 }, we have f3 (x) = k2 − γ2 x2 which points outside of C3 .
– if
k2 γ2
< θ2 + h2 (see similar case shown in Figure 24(b)), when x ∈ {x : −γ1 (θ1 − h1 ) points outside of C3 . When x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, f3 (x) = k2 − γ2 x2 −γ1 x1 . We have that x ∈ {x : x1 > θ1 − h1 , x2 = θ2 − h2 }, f3 (x) = k2 − γ2 (θ2 − h2 ) f3 (x) points inside of C3 .
– If θ2 − h2
θ2 − h2 }, f4 (x) = −γ2 x2 f4 (x) points inside C4 . Then, for every initial condition z(0, 0) ∈ C4 , solutions will reach C1 in finite time. From the above analysis, we have that: 1) from C2 trajectories go to either C1 or C3 ; 2) from C3 trajectories go to either C1 or C4 ; 3) from C4 trajectories go to C1 . Then, trajectories eventually enter C1 , which, using the arguments above, implies that the equilibrium point z2∗ is globally asymptotically stable. 49
Similarly, for case 4 in Table 1, it can be proven that z2∗ is globally asymptotically stable. Since stability of z2∗ was proven in the proof of case 2 in Table 1, now we check the vector fields on the boundaries of C1 when θ1 − h1 < γk11 < θ1 + h1 , θ2 + h2 < γk22 < θ2max (see Figure 27). x2 q1 = 0
q2 = 0
θ2max
C4
C3
C1
C2
θ2 + h2 θ2 − h2
x1 θ1−h1 θ1+h1
θ1max
Figure 27: The vector field at the boundaries of C1 , when θ1 − h1 < k2 < θ2max . γ2
• When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = θ1 − h1
0, we have that f2 (x) points inside of C2 . Then, once a trajectory enters or starts from C2 , it will stay or never leave C2 Then, the set C2 is forward invariant. Since the equilibrium point z1∗ belongs to C2 , every trajectory reaching or starting from C2 converges to z1∗ . For z(0, 0) ∈ C1 , we check the vector fields on the boundaries of C1 (see Figure 31(b)). k1 − γ1 (θ1 + h1 ) . Since • When x ∈ {x : x1 = θ1 + h1 , 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 θ1 − h1 < γk11 < θ1 + h1 , we have that f1 (x) points inside C1 . k1 − γ1 x1 . Since • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = θ2 + h2 }, f1 (x) = −γ2 (θ2 + h2 ) θ1 − h1 < γk11 < θ1 + h1 , we have that f1 (x) points inside C1 . k1 − γ1 (θ1 + h1 ) . f1 (x) is • When x ∈ {x : 0 ≤ x1 ≤ θ1 + h1 , x2 = 0}, f1 (x) = 0 tangent to the boundary of C1 . k1 . Since k1 > 0, we have • When x ∈ {x : x1 = 0, 0 ≤ x2 ≤ θ2 + h2 }, f1 (x) = −γ2 x2 that f1 (x) points inside C1 . 53
x2
x2
q1 = 1
q1 = 0
q2 = 0
q2 = 0
θ2max
θ2max C4
C3
C4
C3
C1
C2
θ2 + h2 θ2 − h2
θ2 + h2 θ2 − h2 C2
C1
x1
x1 θ1−h1 θ1+h1 (a)
θ1−h1 θ1+h1 (b)
θ1max
θ1max
Figure 31: when θ1 −h1 < kγ11 < θ1 +h1 ,0 < γk22 < θ2 +h2 : (a) the vector field at the boundaries of C2 , (b) the vector field at the boundaries of C1 . Then, once a trajectory enters or starts from C1 , it will stay or never leave C1 . Then, the set C1 is forward invariant. Since the equilibrium point z2∗ belongs to C1 , every trajectory reaching or starting from C1 converges to z2∗ . For initial points in C4 , the vector field of the boundary at C4 is shownin Figure 23 (b). −γ1 x1 . Then, we have When x ∈ {x : 0 < x1 < θ1 + h1 , x2 = θ2 − h2 }, f4 (x) = −γ2 (θ2 − h2 ) that f4 (x) points outside C4 and every solution leaves C4 by jumping into C1 . Then, for every initial condition z(0, 0) ∈ C4 , solutions will reach C1 in finite time. As z2∗ ∈ C1 , so the trajectory will stay in C1 and converge to z2∗ . If z(0, 0) ∈ C3 , the vector field at the boundary of C3 is similar as that shown in Figure 24. 1. If the parameters are in the range of θ1 − h1 < kγ11 < θ1 + h1 , 0 < γk22 < θ2 − h2 (see max similar case shown in Figure 24(a)), when x ∈ {x : θ1 − h1 < x1 < θ1 , x2 = θ2 − h2 }, −γ1 x1 . We have that f3 (x) points outside of C3 . When x ∈ {x : f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , which points outside x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, we have f3 (x) = k2 − γ2 x2 C3 . Depending on which jump set the solution hits, two possibilities of the equilibrium points exist. • If the trajectory hits the set x ∈ {x : x1 = θ1 − h1 , x2 > θ2 − h2 }, which leads to an update of q1 , the solution will jump into C4 . For this case, the trajectory will reach C1 in finite time and converge to z2∗ . • If the trajectory hits the set x ∈ {x : x1 > θ1 − h1 , x2 = θ2 − h2 }, q2 is updated to 0 from 1, the trajectory will enter C2 , and converge to the equilibrium point z1∗ , which is inside of C2 . 54
2. If the parameters are in the range of θ1 − h1 < kγ11 < θ1 + h1 , θ2 − h2 < kγ22 < θ2 + h2 (see similar 24(b)), when x ∈ {x : x1 > θ1 + h1 , x2 = θ2 − h2 }, case shown in Figure −γ1 x1 . We have that f3 (x) points inside of C3 . When x ∈ {x : f3 (x) = k2 − γ2 (θ2 − h2 ) −γ1 (θ1 − h1 ) , which points outside x1 = θ1 − h1 , x2 ≥ θ2 − h2 }, we have f3 (x) = k2 − γ2 x2 of C3 . Thus, when the trajectory hits the set {x : x1 = θ1 − h1 , x2 > θ2 − h2 }, a jump will occur, the solution will jump into C4 and finally enter C1 and converge to z2∗ .
55
