Properties of Isostables and Basins of Attraction of Monotone Systems
Abstract
In this paper, we investigate geometric properties of monotone systems by studying their isostables and basins of attraction. Isostables are boundaries of specific forwardinvariant sets defined by the socalled Koopman operator, which provides a linear infinitedimensional description of a nonlinear system. First, we study the spectral properties of the Koopman operator and the associated semigroup in the context of monotone systems. Our results generalize the celebrated PerronFrobenius theorem to the nonlinear case and allow us to derive geometric properties of isostables and basins of attraction. Additionally, we show that under certain conditions we can characterize the bounds on the basins of attraction under parametric uncertainty in the vector field. We discuss computational approaches to estimate isostables and basins of attraction and illustrate the results on two and four state monotone systems.
I Introduction
In many applications, such as economics [1] or biology [2], linear dynamical systems have states, which take only nonnegative values. For example, protein concentrations in biology are always nonnegative. These systems are called positive and have received a considerable attention in the context of systems theory [3], [4], model reduction [5, 6], distributed control [7], [8], etc. Analysis of such systems is facilitated by employing the celebrated PerronFrobenius theorem (cf. [9]), which establishes strong spectral properties of the drift matrix in a linear positive system.
In the nonlinear setting positive systems have a couple of generalizations, with the most known being cooperative monotone systems (cf. [2]). More specifically positive systems generalize to cooperative with respect to the positive orthant systems. Similarly to the linear case, these nonlinear systems generate trajectories (or flows) which for every time are increasing functions in every argument with respect to the initial state. With a slight abuse of notation, we will refer to cooperative monotone systems simply as monotone. In [2], it was briefly mentioned that the flow of a monotone system can be seen as a positive operator. Hence the authors argued that an operator extension of the PerronFrobenius theorem, which is called the KreinRutman theorem [10], can be applied. However, the investigation into spectral properties of these operators lacked due to absence of a welldeveloped theory of spectral elements of such operators. This gap was filled by the development of the socalled Koopman operator and the associated semigroup (cf. [11, 12, 13]).
The Koopman semigroup is a linear infinitedimensional representation of a nonlinear dynamical system. In many practical situations, the spectral elements of this semigroup such as eigenfunctions (which can be seen as infinite dimensional eigenvectors) and eigenvalues can be computed [14, 15]. As in the linear case, the dominant eigenfunctions (i.e., eigenfunctions corresponding to the eigenvalue with the largest real part) offer an insight into the dynamics of the system. In [14], it is shown that the level sets of the dominant eigenfunction (called isostables) contain the initial states of trajectories that converge synchronously toward the equilibrium. The level set at infinity is the boundary of the basin of attraction of the equilibrium. To summarize, the dominant eigenfunctions and eigenvalues contain the information about the asymptotic behavior of the system.
In this paper, we first study the properties of the Koopman semigroup associated with a monotone system. We show that for stable systems monotone with respect to the nonnegative orthant the dominant eigenfunction is an increasing function in every argument for all states in the basin of attraction. Hence we offer yet another version of the PerronFrobenius theorem for nonlinear systems (see [16] for other results on the subject). We use our version of the PerronFrobenius theorem to study geometric concepts of isostables and basins of attraction for monotone systems. Using this statement we provide a straightforward proof of a known result stating that a basin attraction of a monotone system is a union of the socalled orderintervals (cf. [17, 18, 19]). We refine this result by showing that the sublevel sets of the dominant eigenfunctions are unions of orderintervals, as well.
We proceed our theoretical development by considering the basins of attraction of bistable monotone systems under parametric uncertainty in the vector field. We show that if the system is monotone with respect to parameter variations, then it is possible to estimate inner and outer bounds on the basin of attraction of the uncertain system. In this instance, instead of the spectral theory we use monotone control systems theory. Moreover, it appears that the behavior of the eigenfunctions and isostables under parametric uncertainty is more complicated in comparison with the behavior of basins of attraction.
The rest of the paper is organized as follows. In Section II, we cover the main properties of the Koopman operator and monotone systems. In Subsection IIIA we obtain spectral properties of monotone systems, while deriving some properties of the isostables in Subsection IIIB. We investigate the behavior of basins of attraction of monotone systems under parameter variations in Subsection IIIC. We also discuss different algorithms to compute isostables and the boundary of the basins of attraction in Section IV. We conclude the paper with numerical examples in Section V.
Ii Preliminaries
Throughout the paper we consider parameterdependent systems in the following form
(1) 
where , where and for some integers and . We define the flow map , where is a solution to the system (1) with an initial condition and a parameter . We assume that is continuous in on and Lipschitz in on a every compact subset of for every fixed . When it is clear from the context, we will drop the parameterdependence.
Iia Monotonicity
We will study the properties of the system (1) with respect to a partial order induced by cones. A set is called a positive cone if it is closed under addition and multiplication by a nonnegative scalar and if we have for any . A relation is called a partial order if it is reflexive (), transitive (, implies ), and antisymmetric (, implies ). We define a partial order through a cone as follows: if and only if . We will also write if and , and if . We say that the order is standard if , which, with a slight abuse of notation, we denote as without a subscript. We call a set an orderinterval induced by a cone and denote it as . A set is called pconvex if for every , in such that , and every we have that . Systems whose flows preserve a partial order relation are called monotone systems.
Definition 1
The system is monotone with respect to the cones , if for all , and for all and . The system is strongly monotone with respect to the cones , if it is monotone and holds for all provided , , and either or holds.
A certificate for monotonicity with respect to an orthant is called KamkeMüller conditions [20], where some generalizations of this result may also be found.
IiB Koopman Operator
Spectral properties of nonlinear dynamical systems can be described through an operatortheoretic framework that relies on the socalled Koopman operator. The Koopman operator associated with is an operator acting on the functions (also called observables). The Koopman operator generates the semigroup
(2) 
where is the composition of functions and is a solution to the considered system. The Koopman semigroup is linear [11], hence it is natural to study its spectral properties. In particular, the eigenfunctions of the Koopman semigroup are defined as the functions satisfying
(3) 
and is the associated eigenvalue. We can also obtain a very useful expression
(4) 
If is analytic, and if the eigenvalues of the Jacobian of the vector field at the equilibrium are distinct and for all , then the flow of the system can be expressed (at least locally) through the following expansion:
(5)  
where is the space of nonnegative integers, are the right eigenvectors corresponding to , the vectors are the socalled Koopman modes (cf. [11, 15] for more details). In the case of a linear system with matrix having the left eigenvectors , the eigenfunctions are equal to and the expansion (5) has only the finite sum. A similar (but lengthy) expansion can be obtained even if the eigenvalues are not distinct and have linearly dependent eigenvectors (cf. [21]).
Let be such that , , then the eigenfunction , which we call dominant, can be computed using the socalled Laplace average (cf. [14])
(6) 
For all that satisfy and , the Laplace average is equal to up to a multiplication with a scalar. Therefore, without loss of generality, we will only consider socalled increasing functions as observables. A function is called increasing with respect to if for all , such that we have . The eigenfunctions with (nondominant eigenfunctions) are generally harder to compute and are not considered in the present study.
The eigenfunction captures the dominant (i.e. asymptotic) behavior of the system. Moreover for , it follows from (5) that the trajectories starting from the boundary of the set share the same asymptotic evolution
Naturally, these sets are important for understanding the dynamics of the system. The sets are called isostables, and contain the initial conditions of trajectories that converge synchronously toward the equilibrium. In the neighborhood of , the isostables are parallel hyperplanes (if ). We will also employ the notations
for real and nonnegative . A more rigorous definition of isostables and more information can be found in [14].
Iii Geometric Properties of Monotone Systems
Iiia Spectral Properties of Monotone Systems
We start by establishing the properties of the dominant eigenfunctions and eigenvalues of monotone systems.
Proposition 2
Assume that the system admits a stable equilibrium with a basin of attraction . Let be the eigenvalues of the Jacobian at such that for all , and let be corresponding right eigenvectors. Let be the dominant eigenfunction with . Then:
(i) if the system is monotone with respect to on , then is real and negative, for all , satisfying , and . Additionally, we have that for all .
(ii) if the system is strongly monotone with respect to on , then is simple, real and negative, for all , satisfying , and .
Both conditions (i) and (ii) are only necessary and not sufficient, which is consistent with the linear case. The proof of this result is similar to the proof of the main result in [22] concerning spectral properties of the socalled eventually monotone systems. If the vector field is analytic, then the second condition is necessary and sufficient for socalled strong eventual monotonicity.
We can view this result through the prism of positive operator theory. An operator acting on a normed space is called positive (or invariant) with respect to a cone , if . It can be shown that the Koopman semigroup associated with a monotone dynamical system is a positive operator with respect to the cone of increasing functions at every time . Hence the following result may be seen as a nonlinear version of the PerronFrobenius theorem or a version of the KreinRutman theorem for the Koopman operator. We note that a semigroup associated with a nonlinear system was studied from the positive operator viewpoint in the context of diffusion processes in [23].
There were other successful attempts to extend the PerronFrobenius theorem to nonlinear systems (cf. [24, 16, 25]), which however are based on applications of PerronFrobenius arguments to nonlinear maps, while our approach is rather based on their infinitedimensional extension. Moreover, our main contribution is that our version uses eigenfunctions, which give additional insight into the dynamics of the system and can be estimated using databased, algebraic or simulation methods (cf. [15]). Nevertheless, there are similarities, and an investigation into a connection between these results is one of the future work directions.
IiiB Isostables of Monotone Systems
First, we recall a known result on the geometry of basins of attraction and their boundary . This result is briefly mentioned in [2, 18, 17] in the case of bistable monotone systems. The following version requires monotonicity on the basin of attraction only (for a complete proof see e.g. [26]).
Proposition 3
Let the system have an asymptotically stable equilibrium with a domain of attraction and be monotone on with respect to . Then for any , the orderinterval is a subset of .
We proceed by showing similar properties of the sets and isostables .
Proposition 4
Let the system have a stable hyperbolic equilibrium with a domain of attraction and be monotone with respect to on with . Then the following statements hold:
(i) for any , , the orderinterval is a subset of ,
(ii) the manifolds , do not contain , such that ;
(iii) if the system is strongly monotone then the manifold does not contain , such that .
The proof of this result is identical to the case of (strongly) eventually monotone systems [22]. It is important to note that the boundary can contain two points , such that , but in this case these points belong to different manifolds , . The point (ii) can also be shown for , hence the boundaries , do not contain , such that as well. The point (iii) does not generalize in a straightforward manner to the case . This is however shown in [17, 18] under additional assumptions.
IiiC Basins of Attraction of Monotone Systems Subject to Parameter Variations
First we need to discuss the following assumptions.

The system is a bistable system on with two stable equilibria , for all . Assume also that the system is monotone in and .

Parameters take values from a compact set with vectors , such that

The following holds for all :

.
Assumption A1 defines a monotone bistable system, which is also monotone with respect to parameter variations. In order to simplify the presentation we consider only the standard orders for monotonicity in the initial conditions and parameters (which can be verified by means of Proposition 1), but our results hold for other orderings as well. In order to obtain any meaningful bounds on the basins of attraction, we need to assume that the set is compact and ordered, which is done in Assumption A2. We also assume that parameter variations are small enough so that no bifurcation occurs under these parameter variations. Moreover, we assume that the parameter variations do not affect the locations of the equilibria in a significant way, which is done in Assumption A3. This assumption and Assumption A4 are technical and perhaps can be avoided^{1}^{1}1For example, we can modify Assumption A3 to require the sets , to be nonempty, however, this modified assumption is harder to check, but they are straightforward to satisfy and will simplify the presentation of results. The assumptions on the location of stable equilibria and the basins of attraction are schematically depicted in Figure 1. Assumptions A1, A2, and A4 are straightforward to check, while Assumption A3 can only be verified after basins of attraction for and are computed. Using these assumptions we can derive the following result.
Theorem 1
Consider the system under the assumptions A1 – A4. Then and for all . In particular, we have
(7) 
Proof: i) We need to show that for all , we have that . But first we will verify that for all . Due to monotonicity for all we have that
where is equal to for all , while converges to since due to Assumption A3. Hence, for all , such that . Similarly we can show that for all . Together with Assumption A4, this implies that for all , .
ii) Let , due to monotonicity we have that
for , where converges to with . Now, by point i), we conclude that cannot converge to , hence it converges to . This implies that for all .
iii) The inclusion follows from the fact that for all . While (7) is now straightforward to show.
Theorem 1 implies that we can predict a direction of change in the basins of attraction subject to parameter variations if we check a couple of simple conditions. In many applications (e.g. synthetic biology), socalled toggle switches are used, which are bistable systems. The larger the domain of attraction, the more robust is the toggle switch with respect to the corresponding equilibrium. Therefore this result can be valuable for design purposes in applications, where we need to ensure operational robustness of a stable equilibrium in a toggle switch subject to intrinsic and/or exogenous noise.
Iv Tools for Computing Isostables and Basins of Attraction
In control theory, a goto approach to compute forwardinvariant sets (not only basins of attraction) of dynamical systems is sumofsquares (SOS) programming. This approach can be applied to systems with polynomial vector fields and also can provide estimates on the basins of attraction for such systems. There exists a number different SOSbased methods for estimation of basins of attraction, see e.g. [27] and [28] and the references within.
Another option is to compute the function , which provides the isostables and the basin of attraction. In the case of a polynomial vector field, we can formulate the computation of as an infinite dimensional linear algebraic problem using (4). Hence, we can provide an approximation of using linear algebra by parametrizing with a finite number of basis functions (cf. [15]). On another hand, we can estimate directly from data using dynamic mode decomposition methods (cf. [29, 30]). These two options provide extremely cheap estimates of . In fact, the algebraic methods (as demonstrated in [15]) also provide good estimates on basins of attraction. However, we cannot typically compute estimates with an excellent approximation quality, which comes as tradeoff for fast computations.
If the quality of approximation is very important we can simulate a number of trajectories with initial points on a mesh grid of the statespace and compute in these points using Laplace averages (6). After that different interpolation or machine learning methods can be applied to estimate the dominant eigenfunction.
We stress that computing the eigenfunction using Laplace averages is not necessarily more computationally difficult than computing the basin of attraction using sumofsquares. Moreover, the Laplace average method can be applied to systems with nonpolynomial vector fields and is wellsuited to high dimensions. This is due to the fact that the simulation time of a nonlinear system can significantly be reduced by model reduction (timescale separation) methods. In contrast to SOSbased methods, however, the Laplace average method does not generally provide a certificate, that is we cannot guarantee that the estimated set is an inner or an outer approximation of a basin of attraction.
In this paper, we use Laplace averages in our computational algorithm. We use an adaptive grid and thus lower computational complexity by exploiting the facts that the eigenfunction is increasing and only one level set is needed. Finally, monotonicity allows us to guarantee that we compute an inner and an outer approximation of the basin, hence providing a certificate.
Iva Algorithm for Computation of a Level Set of an Increasing Function
First we derive an algorithm for computation of the level sets of an arbitrary increasing with respect to function . Our algorithm can be extended to an arbitrary positive cone in , however, in order to simplify the presentation we focus on the ordering induced by a positive orthant. We compute the level sets within a hypercube . In particular, we compute a set of points , which lies in , and a set of points , which does not lie in . The set (respectively, ) will not contain points , such that . Since is an increasing function this will allow to build a piecewise constant inner approximation (respectively, an outer approximation) of the sublevel set. In order to do so we will use the following oracle
(8) 
We will generate the points on the face of with randomly, while choosing the point greedily. In order to do so efficiently, we need to modify a hypercube in such a way that the curve does not intersect the faces and . Therefore, first we need to increase sufficiently such that all vertices of the face with return the value of the oracle equal to one. At the same time, we need to lower the face with significantly so that all the vertices of this face return equal to zero. Due the constraints on the function we can have and , which is typically the case in the statespace of biological applications. In this case, we need to move the points , and shrink in order to achieve the same goal. We can do so, for example, by computing the level set with and a pick a point with the smallest norm.
Our algorithm exploits the monotonicity property of in the following manner. If a sample is such that , then for all we have . Similarly, if a sample is such that , then for all we have . Therefore, we need to keep track of the largest (in the order) samples with , the set of which we denote , and the smallest (in the order) samples with , the set of which we denote . Hence, in order to approximate the function at every step we generate new samples such that . After generating the samples we compute the oracle values and add new samples to the set if , or to the set if . After new samples have been added there might exist , in (respectively, ) such that (respectively, ). In these cases all such samples must be removed from (respectively, ). We call this procedure pruning of the sets.
All that is left is to explain the procedure of generation of new samples. First of all, let with the Lebesgue measure . The measure of the whole space, which is the defined above hypercube , is its volume . The value represents the current error of the algorithm. In order to sample from the actual level set the measure has to be equal to zero.
We consider two ways of generating samples: random and greedy. In “random” generation we randomly generate points in the set . Let be a sample, which we need to generate. First, we generate the scalars , , on the edge the hypercube with . We do it by using a probability distribution with the support on the whole edge. Then we compute the limits for the generation of such that . Next, we generate using the same distribution , while adjusting its support. We generate samples in this manner.
The “greedy” generation is meant to exploit more information about and decrease the value , which represents the error of the algorithm. For every sample in we find points from such that and such that there is no other satisfying . Then we form hypercubes with vertices and and choose with the largest volume. We denote these hypercubes . Having done it for every sample in we compute hypercubes with largest volumes. After that we generate samples in each of the hypercubes using a probability distribution . If no information about the system is known then we choose , as uniform distributions. Otherwise, we can choose a Beta distribution so that the median approximates the most probable location of the curve.
The procedure is summarized in Algorithm 1. Note that our sample generation guarantees that the measure is nonincreasing with generation of new samples, while random sampling guarantees that the value eventually converges to zero with the rate of convergence depending on specific functions . However, in our examples we observed exponential empirical convergence, which can be explained by the fact that our algorithm can be seen as a bisection procedure using random sampling.
IvA1 Oracles for Computation of Isostables and Basins of Attraction
Without loss of generality we assume that the stable equilibrium lies at the origin. Since the sets and are level sets of some increasing function according to Proposition 4, we can use Algorithm 1 to compute them. In order to do so, we can use the following oracle:
(9) 
where the value of the eigenfunction is computed using Laplace averages (6), and is a large enough time. We need the second condition to make sure that the points also lie in . Even though it is theoretically unlikely to have a point with a finite (nonzero) not lying in , such situations can occur numerically. If we need to compute , then we can drop the first condition , since in this case .
We can also apply this algorithm to the computation of the socalled switching separatrix defined in the context of pulsecontrolled bistable systems [26], since it is also the level set of an increasing function.
V Numerical Examples
Va A TwoState LacITetR Switch
The genetic toggle switch is composed of two mutually repressive genes LacI and TetR and was a pioneering genetic system for synthetic biology [31]. We consider the following model of the toggle switch
(10) 
where all . The states represent the concentration of proteins LacI and TetR, whose mutual repression is modeled via a rational function. The parameters and model the basal synthesis rate of each protein. The parameters and are degradation rate constants, and , describe the strength of mutual repression. The parameters , are called Hill coefficients. By means of Proposition 1 we can readily check that the model is monotone on for all nonnegative parameter values with respect to the orthant . Moreover, the model is monotone with respect to all parameters but , . The stable equilibrium has the state “switched on” (i.e., is much larger than ), while has the state “switched on” (i.e., is much larger than ).
In order to evaluate Theorem 1 we consider the set of admissible parameters , where
We compute the isostables with , , for systems with parameters , , , where is the boundary of the basin of attraction , and the matrix of parameters is equal to:
The computational results are depicted in Figure 2. These results suggest that for all parameter values the manifold will lie between the manifolds and . It appears that also lies between the manifolds and , however, in a different order. This change of order and continuity of implies that there exists an such that at least two manifolds , , intersect. This case is also depicted with . This observation implies that is not an increasing function in . This is consistent with the linear case, where changes in the drift matrix , will simply rotate the hyperplane , where is the left dominant eigenvector of .
VB A Four State Toggle Switch.
We can compute the basins of attraction and isostables for dimensional models, as well. Naturally, we can visualize only D crosssections of basins of attraction of dimensional systems. We consider again a toggle switch, while modeling also mRNA concentrations activating their corresponding proteins. The resulting model is as follows
(11) 
where the parameters have the following values
Again, by means of Proposition 1, we can check that the model is monotone on for all nonnegative parameter values with respect to the orthant . We will plot a crosssection with . In Figure 3, we plot isostables , , and . We also plot the boundary of the domain of attraction.
It is noticeable that the sets do not depend on the state in a significant manner. Hence the state can be reduced using timescale separation methods. Similarly, we can reduce . This is in an agreement with “a rule of thumb” in proteinmRNA interactions, which states that mRNA dynamics are much faster than protein dynamics and hence can be reduced. Moreover, for every fixed , the shapes of isostables are similar to the shape of the isostables of a two state toggle switch.
Vi Conclusion
In this paper, we study geometric properties of monotone systems such as properties of isostables and basins of attraction. Isostables contain the initial states of trajectories that converge synchronously toward the equilibrium and can be computed using the Koopman operator framework. Isostables are boundaries of forwardinvariant sets of dynamical systems, and hence serve as convenient refinement of the basin of attraction notion. We show that the isostables and the boundary of basin of attraction of monotone systems share the same geometric properties. For example, these manifolds do not contain strictly comparable points in the order. Our derivations are based on a positivity result for the Koopman semigroup, which is related to the statement of the PerronFrobenius theorem. We then focus on properties of basins of attraction of bistable systems. We show how to estimate basins of attraction of monotone parameterdependent systems subject to parametric uncertainty, and illustrate our results on numerical examples. Future work includes more practical implications of our results in synthetic biology, and a further study of the properties of our algorithm such as rate of convergence, efficiency and scalability.
References
 [1] W. W. Leontief, Inputoutput economics. Oxford University Press on Demand, 1986.
 [2] M. Hirsch, H. Smith, et al., “Monotone dynamical systems,” Handbook of differential equations: ordinary differential equations, vol. 2, pp. 239–357, 2005.
 [3] P. G. Coxson and H. Shapiro, “Positive input reachability and controllability of positive systems,” Linear Algebra and its Applications, vol. 94, pp. 35–53, 1987.
 [4] L. Farina and S. Rinaldi, Positive linear systems: theory and applications. John Wiley & Sons, 2010.
 [5] A. Sootla and A. Rantzer, “Scalable positivity preserving model reduction using linear energy functions,” in Proc IEEE Conf Decision Control, Dec. 2012, pp. 4285–4290.
 [6] C. Grussler and T. Damm, “A symmetry approach for balanced truncation of positive linear systems,” in 51st IEEE Conference on Decision and Control, Maui, Hi, USA, Dec. 2012, pp. 4308–4313.
 [7] A. Rantzer, “Scalable control of positive systems,” European Journal of Control, vol. 24, pp. 72–80, 2015.
 [8] T. Tanaka and C. Langbort, “The bounded real lemma for internally positive systems and hinfinity structured static state feedback,” IEEE transactions on automatic control, vol. 56, no. 9, pp. 2218–2223, 2011.
 [9] A. Berman and R. J. Plemmons, Nonnegative Matrices in the Mathematical Sciences. SIAM, 1994, vol. 9.
 [10] M. G. Krein and M. A. Rutman, “Linear operators leaving invariant a cone in a Banach space,” Amer. Math. Soc. Translation, p. 128, 1950.
 [11] I. Mezić, “Spectral properties of dynamical systems, model reduction and decompositions,” Nonlinear Dynam, vol. 41, no. 13, pp. 309–325, 2005.
 [12] M. Budišić, R. Mohr, and I. Mezić, “Applied koopmanism,” Chaos, vol. 22, no. 4, p. 047510, 2012.
 [13] I. Mezic, “Analysis of fluid flows via spectral properties of the Koopman operator,” Annual Review of Fluid Mechanics, vol. 45, pp. 357–378, 2013.
 [14] A. Mauroy, I. Mezić, and J. Moehlis, “Isostables, isochrons, and Koopman spectrum for the action–angle representation of stable fixed point dynamics,” Physica D, vol. 261, pp. 19–30, 2013.
 [15] A. Mauroy and I. Mezic, “Global stability analysis using the eigenfunctions of the Koopman operator,” IEEE Tran Autom Control, 2016, In press.
 [16] B. Lemmens and R. Nussbaum, Nonlinear PerronFrobenius Theory. Cambridge University Press, 2012, vol. 189.
 [17] H. Smith and H. Thieme, “Stable coexistence and bistability for competitive systems on ordered banach spaces,” Journal of Differential Equations, vol. 176, no. 1, pp. 195–222, 2001.
 [18] J. Jiang, X. Liang, and X.Q. Zhao, “Saddlepoint behavior for monotone semiflows and reaction–diffusion models,” J. Differential Equations, vol. 203, no. 2, pp. 313–330, 2004.
 [19] B. S. Rüffer, H. Ito, and P. Dower, “Computing asymptotic gains of largescale interconnections,” in IEEE Proc Conf Decision Control, 2010, pp. 7413–7418.
 [20] D. Angeli and E. Sontag, “Monotone control systems,” IEEE Trans Autom Control, vol. 48, no. 10, pp. 1684–1698, 2003.
 [21] I. Mezic, “On applications of the spectral theory of the koopman operator in dynamical systems and control theory,” in IEEE Conf Decision Control, 2015, pp. 7034–7041.
 [22] A. Sootla and A. Mauroy, “Operatortheoretic characterization of eventually monotone systems,” Provisionally accepted for publication in Automatica, March 2016, http://arxiv.org/abs/1510.01149v2.
 [23] I. Herbst and L. Pitt, “Diffusion equation techniques in stochastic monotonicity and positive correlations,” Probab Theory Rel, vol. 87, no. 3, pp. 275–312, 1991.
 [24] S. Gaubert and J. Gunawardena, “The PerronFrobenius theorem for homogeneous, monotone functions,” Transactions of the American Mathematical Society, vol. 356, no. 12, pp. 4931–4950, 2004.
 [25] D. Aeyels and P. De Leenheer, “Extension of the Perron–Frobenius theorem to homogeneous systems,” SIAM Journal on Control and Optimization, vol. 41, no. 2, pp. 563–582, 2002.
 [26] A. Sootla, D. Oyarzún, D. Angeli, and G.B. Stan, “Shaping pulses to control bistable systems: Analysis, computation and counterexamples,” Automatica, vol. 63, pp. 254–264, Jan. 2016.
 [27] D. Henrion and M. Korda, “Convex computation of the region of attraction of polynomial control systems,” IEEE Tran Autom Control,, vol. 59, no. 2, pp. 297–312, 2014.
 [28] G. Valmorbida and J. Anderson, “Region of attraction analysis via invariant sets,” in Proc Am Control Conf, 2014, pp. 3591–3596.
 [29] P. J. Schmid, “Dynamic mode decomposition of numerical and experimental data,” Journal of Fluid Mechanics, vol. 656, pp. 5–28, 2010.
 [30] J. H. Tu, C. W. Rowley, D. M. Luchtenburg, S. L. Brunton, and J. N. Kutz, “On dynamic mode decomposition: Theory and applications,” J Comput Dynamics, vol. 1, no. 2, pp. 391 – 421, December 2014.
 [31] T. Gardner, C. R. Cantor, and J. J. Collins, “Construction of a genetic toggle switch in escherichia coli,” Nature, vol. 403, pp. 339–342, 2000.