On convergence for hybrid models of gene regulatory networks under polytopic uncertainties: a Lyapunov approach

Hybrid models of genetic regulatory networks allow for a simpler analysis with respect to fully detailed quantitative models, still maintaining the main dynamical features of interest. In this paper we consider a piecewise affine model of a genetic regulatory network, in which the parameters describing the production function are affected by polytopic uncertainties. In the first part of the paper, after recalling how the problem of finding a Lyapunov function is solved in the nominal case, we present the considered polytopic uncertain system and then, after describing how to deal with sliding mode solutions, we prove a result of existence of a parameter dependent Lyapunov function subject to the solution of a feasibility linear matrix inequalities problem. In the second part of the paper, based on the previously described Lyapunov function, we are able to determine a set of domains where the system is guaranteed to converge, with the exception of a zero measure set of times, independently from the uncertainty realization. Finally a three nodes network example shows the validity of the results.


Introduction
In the last few years control theory tools have been extensively used in biology, to both understand natural biological systems or design new ones to perform specific tasks (Blanchini et al. 2018;Qian et al. 2018). Within a cell, in fact, multiple regulatory mechanisms coexist (Alon 2007;Freeman 2014) and among these, transcriptional regulation involving genes and transcription factors plays a crucial role. Transcription factors are proteins which can either activate or inhibit the transcription of different genes, which in turn can produce other transcription factors, determining a set of relations described by a gene regulatory network (GRN) (Alon 2007). Different modelling approaches are used to mathematically describe a GRN, depending on whether the analysis to perform is either quantitative or qualitative (Karlebach and Shamir 2008). Continuous models, derived from the study of chemical reaction networks and quasi-steady state approximation of faster dynamics, usually involve Hill functions to describe regulatory interactions and are particularly well suited for a quantitative analysis (Alon 2007;Le Novère and Nicolas 2015;Murray and Del Vecchio 2014). However, a quantitative approach may not always be possible as biological systems are inherently uncertain and the measurement of key quantities may be affected by noise or impossible to take; for these reasons, qualitative approaches have been developed in literature. Asynchronous boolean networks have been studied to analyse complex biological systems (see for example Tournier and Chaves 2009 and references therein), as many interactions and system quantities can be described by logic variables (e.g. HIGH or LOW concentration, protein production is ACTIVE or INACTIVE, etc.). Despite being a valid and simple method, completely discrete analysis can lose track of some dynamical behavior (Saadatpour and Albert 2016) and it may be necessary to consider hybrid models. The hybrid modelling approach was introduced in Glass and Kauffman (1973) and has been studied and adapted by many authors (see for example Casey et al. 2006;De Jong et al. 2003, 2004Ropers et al. 2006;Cummins et al. 2016;Gedeon 2020). Based on the step approximation of steep Hill functions (Alon 2007), this approach gives rise to a piecewise affine (PWA) model. In Casey et al. (2006) the authors studied this model and gave stability conditions of equilibria based on a State Transition Graph (STG), a graph that qualitatively describes the system trajectories. The STG is unchanged for a large range of parameters, making this kind of analysis attractive when the system's knowledge is partial. Despite being a valid tool, there are a few limitations in the use of STGs. These graphs in fact, besides introducing spurious qualitative trajectories (i.e. paths in the graph that do not correspond to trajectories of the original system) (Casey et al. 2006), do not give clear answers when cycles are present (Grognard et al. 2007), making it difficult to distinguish between damped oscillations and limit cycles. This motivates the development of tools that can fill these gaps, by complementing these qualitative tools, exploit-ing additional quantitative information from the system. One of the most common approach in this sense is to consider Lyapunov functions, which can help to discriminate between different kind of dynamical behaviors, otherwise indistinguishable from the graph, and the use of which is even suggested in the Conclusion of Casey et al. (2006), and which are the object of this work. The use of Lyapunov functions in the study of biochemical networks is well known in literature. In Ali Al-Radhawi and Angeli (2016) the authors studied chemical reaction networks stability and control, using piecewise linear in rates Lyapunov functions, in Blanchini and Giordano (2014) the authors employed piecewise linear Lyapunov functions to assess structural properties of biological networks, while in Chesi and Hung (2008) the authors, based on a linear matrix inequalities (LMIs) framework, used polynomial Lyapunov functions to study the stability of genetic regulatory networks with SUM regulatory functions. In our previous work (Pasquini and Angel 2019) we considered the aforementioned PWA model of a GRN and, using the model structure and the information from the STG, we developed an LMI framework to find a Lyapunov function for the system, in order to assess its convergence properties, even in the presence of cycles in the STG. In that context we assumed complete knowledge of the system's parameters, an assumption that in general does not hold. However, Lyapunov methods are used also in the uncertain setting and a particular case is that of polytopic uncertainties for linear systems, in which the system's matrix is the convex combination of a set of known matrices, with the weights of this combination being unknown. Two approaches are generally employed in literature. The first one consists in searching for a Lyapunov function which is common between all the systems obtained by considering the realizations associated with the vertices of the uncertainty polytope (Liberzon 2003;Lin and Antsaklis 2005). This approach, although allowing to conclude certain stability and stabilizability properties of the system, is challenging and can lead to a conservative solution, as it does not consider how the extremal behaviors are combined in any given uncertainty realization. The other approach, which we will consider in this work to deal with polytopic uncertainties, is to search for a parameter dependent Lyapunov function (PD-LF). In Gahinet et al. (1996) the authors first showed how to convexify the problem of finding an affinely parameter dependent Lyapunov function, through the solution of a set of LMIs. Following this many authors proposed different LMIs framework to deal with linear parameter varying systems stability (Chesi et al. 2004;Neto 1999;Oliveira and Peres 2006) and stabilizability (Lin and Antsaklis 2007;Zhai et al. 2003) through the use of PD-LFs, which are in general polynomially dependent on the uncertain parameters. Although the literature on parameter dependent lyapunov functions for polytopic uncertain systems is vast and still developing, many aspects are not usually considered. In most of the cases where both polytopic uncertainties and switching systems are considered, the problem of stabilizability is addressed (i.e. the choice of a switching signals to stabilise the system). However, in a PWA model of a GRN, the switching is state dependent and cannot be chosen. Moreover multistability, which is a common property of biological systems, is usually not considered; instead global asymptotic stability of a single equilibrium is a common theme in many of these works. For example in Arcak and Sontag (2008) the authors propose techniques based on dissipativity of individual gene subsystems in order to build separable Lyapunov functions that can be used to assess global asymptotic stability of the network. LMI conditions are used to infer the parameter ranges where such conditions can be fulfilled. In contrast, our paper does not focus on global asymptotic stability properties and exploits the structure of PWA models in order to build robust Lyapunov function for global convergence analysis.
In the present paper we consider polytopic uncertainties in piecewise affine models of GRNs, to address the issue that the system's parameters may not be perfectly known (as in qualitative methods), but bounded within preassigned ranges. We provide an LMI framework whose solution describes a parameter dependent piecewise quadratic Lyapunov function (PD-PWQ-LF) for the system, i.e. a Lyapunov function that depends explicitly on the unknown parameters describing the particular uncertainty realization, which will allow us to describe a convergence set for the system, robust with respect to the uncertainty. Our analysis can deal with multistable systems, as information on equilibria location is not explicitly taken into consideration when defining the LMI framework. Moreover it is remarked that even if in literature there are results on parameter dependent Lyapunov functions with more complex and general structures (see for example Chesi et al. 2005), in this work we consider parameter dependent Lyapunov functions which are only affinely dependent on the system's uncertain parameters, as these allow to deal in a straightforward manner with sliding mode monotonicity and continuity on the boundary of regulatory domains, as it will be clear in the following. The paper is structured as follows: in Sect. 2 some mathematical notation and preliminaries are defined, while in Sect. 3 the piecewise affine model and its differential inclusion extension are described, together with a brief recap of our work (Pasquini and Angel 2019) on piecewise quadratic Lyapunov functions for the nominal PWA model. In Sect. 4 the main contribution of the paper is presented: first we describe the considered polytopic uncertain model and the desired form of the Lyapunov function, then, after a discussion on the constraints that need to be enforced by this function, we give the conditions for its existence in the form of an LMI framework (Theorem 1). Finally we prove that, for any realization of the system uncertainty, the set where the system converges is contained in a computable set depending only on the vertices value of the uncertainty polytope. In Sect. 5 we give a numerical example showing the applicability and the validity of the results for a three nodes GRN, while Sect. 6 concludes the paper and give some possible directions for future developments.

Mathematical background
With v ≥ 0 we mean that all the components of v are non-negative. The set of all v ∈ R n such that v ≥ 0 is denoted with R n + . Let W = {w 1 , . . . , w m } be a set of vectors. We define the conic hull of W as cone(W ) A set-valued map H : X → 2 Y is a map that associates a point in X to a subset of Y .
With v∈W we denote the sum over the elements of the countable set W , while with v ∈W we denote the union over the elements of the set W .
Given a matrix M ∈ R n×n , M i j denotes the element in its i-th row and j-th column. Let M ∈ R n×n be a symmetric matrix. We denote with M 0(M 0) the fact that M is positive definite (semidefinite), and with M ≺ 0(M 0) the fact that M is negative definite (semidefinite).
A polyhedron P in R n , is the set described by: Equation (1) is called the H -representation of P. Every polyhedron's H -representation can be converted to a V -representation (Avis et al. 2002;Herceg et al. 2013;Iervolino et al. 2017b), namely: where V = v 1 . . . v ν denotes the set of vertices of P, and R = r 1 . . . r ρ the set of its rays, with the notation in (2) meaning that any element in P can be expressed as the sum of an element belonging to conv{V } and an element belonging to cone{R}. Any polyhedron P in R n can be embedded in an higher dimensional cone, called the homogenization cone, i.e. the cone: The homogenization cone of P has the property that P can be obtained as the intersection of the cone with an hyperplane H in R n+1 , thus allowing to give sufficient conditions which guarantee a chosen quadratic function to be either positive or negative definite (semidefinite), inside the original polyhedron P.
In particular we recall the following property from Iervolino et al. (2017a), which we also used in Pasquini and Angel (2019), to express sign definiteness conditions of quadratic functions inside a polyhedron P.
Consider the quadratic function: where: x Let P be a polyhedron andĈ P its homogenization cone, defined above. Then it holds: where Γ is the matrix whose columns are the rays of the coneĈ P , N is any entrywise non-negative matrix and M Ĉ P 0 means that the quadratic function (3) is nonpositive for any x ∈Ĉ P , and consequently in the polyhedron P. For further details the interested reader is referred to Iervolino et al. (2017a).

Hybrid model
In this section the hybrid model considered to describe the GRN dynamics is presented. This section is mainly recalled from (Casey et al. 2006;De Jong et al. 2004;Pasquini and Angeli 2018), to which the reader is referred for a more comprehensive discussion Let C ∈ R n×n be a diagonal matrix with positive entries and let f : R n + → R n + be a piecewise constant function, defined on a box partition of the positive orthant (partition that will be characterised below). We consider the following hybrid model: where x = x 1 . . . x n T ∈ R n + represents the protein concentration vector (with x i being the concentration of the protein P i ), C = diag(c 1 , . . . , c n ), with c i > 0, represents the degradation rates matrix and f (x) = f 1 (x) . . . f n (x) T represents the production rate function, sometimes referred to as the regulation function. Let any axis X i of the positive orthant, be partitioned as: The θ s in (6) are called thresholds. Such partition divides the state space in open boxes, inside which we consider the regulation function f (x) to be constant. These boxes are called regulatory domains, as opposed to the ones called switching domains, in which at least one of the state variables takes on a threshold value. We define the order of a switching domain D s as the number of switching variables in D s i.e. state variables that are on one of their threshold values. The set of regulatory domains is denoted by D R , while the set of switching domains is denoted by D S . Inside each regulatory domain D the production function is considered constant and indicated as f D , hence the dynamics is affine, with the property that any trajectory monotonically converges towards the focal point φ(D) = C −1 f D and, if this is not in cl(D), then the trajectory will eventually leave the domain (Casey et al. 2006).

Remark 1
The above definition of the production function for system (5), is a consequence of the ON-OFF approximation of the gene input functions (Alon 2007). In fact in works like Casey et al. (2006) and De Jong et al. (2004), and our previous work Pasquini and Angeli (2018), the production function is considered as a combination of step functions, the thresholds of which give rise to the same box partition described above defining, de facto, a piecewise constant function. In Plahte et al. (1998) the authors gave a general framework which allows to describe any logical function with the use of steps or sigmoids functions.
The switching domains are the zero measure sets where the production rate function is not uniquely defined and so a further construction due to Filippov (1988) is used. In particular the system is extended to a differential inclusion: for which: in which R(D) denotes the set of regulatory domains adjacent to the switching domain D. A solution (in the sense of Filippov) of system (5) is an absolutely continuous function x(·), satisfying the differential inclusion (7), for almost every t, with the set-valued map F defined as in (8). It is possible that a solution x(·) of the differential inclusion (7), lays for a certain amount of time on a switching domain (i.e. on the surface of discontinuity of f (x)), leading to a so called sliding mode solution.
As explained in Casey et al. (2006) and De Jong et al. (2003, 2004, it is possible to define a state transition graph of the PWA system, a graph that qualitatively characterise the system trajectories in relation to the domains (both regulatory and switching). This construct will appear in a later assumption (Assumption 1), but for more details on how the STG can be constructed and used, the interested reader is referred to the aforementioned literature.
We now introduce a simple example that will help clarify the above notation.

Example: toggle switch
Consider the following piecewise affine model of a toggle switch (the GRN of which is depicted in Fig. 1):  (9) where s − (x i , θ j ) is the step function: and b 10 = 0.02, b 20 = 0.08, b 11 = 3, b 21 = 3.5, c 1 = 0.7, c 2 = 1.1, θ 1 = θ 2 = 2. The toggle switch is one of the simplest gene regulatory networks to consider [and one of the first networks that has been built synthetically (Collins et al. 2000)], in which two proteins P 1 and P 2 act as each other inhibition transcription factor. For certain sets of parameters, the system is known to show bistability, as the two stable configurations are when one of the two proteins concentration is high, while the other one is low.
As can be seen in Fig. 2, the thresholds θ 1 and θ 2 give rise to four regulatory domains (i.e. D 1 , D 2 , D 3 , D 4 ), and five switching domains: four of which of order 1 (e.g. the domain D = ∂ D 1 ∩ ∂ D 2 ) and one of order 2 (i.e. The procedure described in De Jong et al. (2004) generates the state transition graph in Fig. 3. In this graph: the blue nodes are associated to regulatory domains, while the red ones are associated to switching domains (one circle: domains of order 1, two circles: domains of order 2). It is clear from it that the system exhibits bistability, property that is also confirmed by the trajectories in Fig. 4.
The intersection between the closure of the four regulatory domains (i.e. the point (θ 1 , θ 2 )) is a singular equilibrium-in the sense that 0 ∈ F(θ 1 , θ 2 ), with F(x) being the set-valued map defined in (8)

Piecewise quadratic Lyapunov function
Lyapunov functions are a tool that is extensively used in control theory and the study of dynamical systems, allowing the study of stability and convergence properties (Khalil 2002). In Pasquini and Angeli (2018) we presented an LMI framework to find, if one exists, a Piecewise Quadratic Lyapunov function (PWQ-LF) for system (5), which can be formally proved to be eventually non-increasing along any system trajectory (Pasquini and Angel 2019). In particular the Lyapunov function V is defined as: and it is asked to satisfy the following constraints: where T x D is the tangent cone to D in x, ∇V D represents the gradient of the function VD, withD being any regulatory domain adjacent to D, and C D denotes the set of switching domains on which we asked the function V to be continuous, namely the ones for which the set of sliding mode directions is non-empty and the ones associated to cycles in the State Transition Graph.
To summarise: constraints (12) and (13) are relative to the monotonicity of the Lyapunov function in regulatory domains and switching domains with sliding modes respectively, while (14) refers to the continuity of the Lyapunov function in particular switching domains-the set of which is referred as C D.
These constraints are enforced through a set of LMIs, which is omitted here for the sake of space and ease of explanation, however, for further details on the description and derivation of this set, the interested reader is referred to our previous works (Pasquini and Angel 2019;Pasquini and Angeli 2018). In particular it can be proven that, for a function V satisfying (12)-(14), the following property holds: Proposition 1 Consider the system (5) and let V : R n + → R be a Lyapunov function satisfying (12)-(14). Then: where μ(S) denotes the Lebesgue measure of the set S.
In the following we refer to property (15) as a result of convergence to zero of dV dt (x(t)) in the sense of measure.
In words, condition (15) states that, given an arbitrarily small ε > 0 and a set Ω where the solution x(·) is contained [e.g. Ω = R n + or Ω = B, with B being a positively invariant set with respect to (5)], then as t → ∞ the solution will be almost surely (in the sense of Lebesgue measure) in the set Ω\{x ∈ Ω |V (x) < −ε}. This means that, by taking values of ε converging to 0, we can infer information on the convergence set for the system. (2019) a few more technicalities are introduced to avoid the convergence to a set where the Lyapunov function is not defined. This has been done exploiting a construction-called natural extension of the Lyapunov function-which is equal to the original function V almost everywhere. More details can be found in the aforementioned paper.

Remark 2 In Pasquini and Angel
The existence of a feasible solution to this set of LMIs, and consequently of a Lyapunov function for the system, is strictly dependent on the system parameters, which are often highly uncertain due to the nature of biological systems. In the next Section the case of uncertainties on the production rate function is considered.

Main contribution
We now introduce polytopic uncertainties on the production rate function, to model the fact that its exact value is unknown and possibly subject to variability.
Let C ∈ R n×n be a diagonal matrix with positive diagonal entries and let f 1 (x), . . . , f L (x) be L piecewise constant production functions, as defined in Sect. 3. Consider the system Σ k , defined as: in which: and S L is the standard simplex of dimension L. In the following we will use the shortened notation: as the dependence of the set from the considered functions f 1 , . . . , f L is implicit. From now on the following Assumption is considered to be satisfied: Assumption 1 All extremal systems Σ 1 , . . . , Σ L have the same State Transition Graph (STG) and the same thresholds.
Remark 3 The first part of Assumption 1 is not too restrictive, in the sense that the STG is generally unchanged under a large range of parameters (see Casey et al. 2006;De Jong et al. 2004), and assuming that the STG is the same for all extremal systems means that these are "qualitatively" similar, which is the starting point of this analysis.
If multiple graphs arise one could, in principle, attempt a partition of the uncertainty set, so as to recover Assumption 1 on each of its element.
On the other hand to consider that thresholds are unchanged among all extremal systems can be artificial. However, considering uncertain thresholds, entails a structural change in the positive orthant partition, with remarkable difficulties in imposing continuity constraints-which will potentially convert in non-linear constraints, given how condition (14) is ultimately enforced (see (Pasquini and Angel 2019, Eq. (43)) or Eq. (42)-(44) of this work). We recognize the need to address this matter in potential future research.

Example: toggle switch with polytopic uncertainties
To better clarify the new notation introduced above, consider again the toggle switch example of Sect. 3.2, but this time allow We can define the four extremal systems: It is easy to verify (see Casey et al. 2006 for example) that all extremal systems share the same STG. The set U 4 1 , of possible system realizations given the uncertainty, can be rewritten in a parametrized form as: with λ ∈ S 4 , or equivalently: By varying λ ∈ S 4 (or η ∈ S 1 × S 1 ) we obtain different realizations σ λ ∈ U 4 1 , with the two stable equilibria moving depending on such parameters, while saddle point (θ 1 , θ 2 ) remains unchanged (following the structure of the system and Assumption 1). In Fig. 5, trajectories for different realizations σ λ ∈ U 4 1 are shown, starting from the same initial points.
In the following we will prove that we can define a set of linear matrix inequalities (LMIs) that, if satisfied, guarantees the existence of a Parameter Dependent Piecewise where V k (x) is a particular PWQ-LF for the k-th extremal system Σ k , will be called an extremal Lyapunov function and will assume the form (11). Formally speaking, given a certain λ ∈ S L , V λ (x) is a PWQ-LF for the system σ λ ∈ U 1 L , and this will hold for all the λ ∈ S L .
The set of LMIs to satisfy will take into account: the conditions to make V k (x) an extremal Lyapunov function for Σ k and the conditions to guarantee that V λ (x) is a PWQ-LF for σ λ .
Remark 4 Asking for V k (x) to be an extremal Lyapunov function for the system Σ k , corresponds to asking (12)-(14) for V k , with respect to Σ k , with the only exception of the constraints on monotonicity along sliding modes (13), that should be substituted as explained in the next section. This is due to the fact that asking for V k to be nonincreasing along any sliding mode solution of Σ k , is not enough to guarantee that V λ in (19) will be non-increasing along any sliding mode solutions of σ λ .

Sliding mode directions description
Consider a switching domain D s . For any system in U L 1 , any potential sliding mode direction f on D s , can be expressed as: where γ ∈ S L·q and q := |R(D s )|. The γ s in (20) are the weights of a convex combination that takes into account all extremal systems, in order to generate candidate sliding mode directions at x ∈ D s . Let I D s be the set of switching variables in D s (i.e. the variables that are on their thresholds in D s ).
The following polyhedron P D s can be defined: wherec is a vector, the i-th component of which is: where c i is the degradation rate of the protein P i and θ i,k is the threshold value assumed by x i in the switching domain D s [see axis partition (6)], while F k is the matrix: in which ( f k D ) I Ds is the vector obtained by selecting only the components from ( f k D ), indexed by I D s , and {D i 1 , . . . , D i q } is the set of regulatory domains adjacent to D s . P D s contains all the γ s that give rise to sliding mode directions and, by construction, it is the same among all the extremal systems and all the systems in U L 1 . Being P D s a subset of the standard simplex S L·q , it is bounded and so can be written as: Following an approach similar to the one we used in Pasquini and Angeli (2018), if every extremal Lyapunov function V k is monotone along the directions obtained by selecting the γ s in P D s , then the function V λ will be monotone along any sliding modes solution.
The above property can be satisfied, by enforcing the following set of LMIs for any extremal LF V k : where Γ D s is the ray matrix of the homogenization coneĈ D s of D s , M k D s is any entrywise non-negative and symmetric matrix and: where P k D and d k D are the matrices P D and d D in (11), for the extremal LF V k , w j is the j-th vertex of P D s and F is the matrix defined as: with: where {D i 1 , . . . , D i q } is the set of regulatory domains adjacent to D s .

Remark 5
The above approach is conservative, in the sense that not every γ ∈ P D s generate an actual sliding mode direction for the systems in U L 1 , but P D s surely contains all of them.

Existence of a PD-PWQ-LF
The following Theorem gives the conditions to construct a PD-PWQ-LF of the form (19).
Theorem 1 Let Σ 1 , . . . , Σ L be L extremal systems as defined in (16). Let V 1 , . . . , V L be their extremal LFs, each one satisfying the LMIs guaranteeing (12) and (14) and the set of LMIs (25). Let, for any regulatory domain D ∈ D R : in which M k j is any non-negative entrywise symmetric matrix, Γ D is the ray matrix of the homogenization coneĈ D and: where: Let σ λ ∈ U L 1 . Then: is a PWQ-LF for σ λ .
Proof Fix a λ ∈ S L . The system σ λ is a PWA system and for V λ to be a PWQ-LF for σ λ , it needs to satisfy the constraints (12), (13) and (14). V λ is a convex combination of the extremal Lyapunov functions and, given Assumption 1, (14) is naturally satisfied for V λ .
Moreover, given (25), every extremal Lyapunov function is monotone along every sliding direction of σ λ , for any switching domain D s , because of the construction of P D s , and this guarantee that V λ satisfies (13) as well. We now only need to prove that (12) is satisfied for V λ with respect to σ λ . Let D be a regulatory domain of (5).
in which: We can rewrite (33) as: Using (34) and the definition of δ f k j D , (35) can be written as: It is easy to show that: with δP k j D defined by (30), and that: Given that λ ∈ S L , it holds: and so using (37) and (38), Eq. (35) can be written as: which, given the constraints (29), immediately gives the following: At this point using the fact that any extremal Lyapunov function V k satisfies (12), it follows that also V λ satisfies (12), proving that V λ is a PWQ-LF for σ λ . The same reasoning is valid for any λ ∈ S L , and this completes the proof.

Remark 6
We refer to the set of LMIs (29) as crossed conditions, as each one of these connects different extremal Lyapunov functions and different extremal systems, inside the same regulatory domain. The need for these conditions is made explicit in equation (39), as the derivative in time of the Lyapunov function V λ D is not only dependent on the derivatives in time of the extremal Lyapunov functions, but also on a set of crossed terms, which are guaranteed to be non-positive if the crossed conditions are satisfied.
These terms are a consequence of the fact that both the system σ λ and the chosen PD-LF V λ are convex combinations of their extremal counterparts.

Extended feasibility problem definition
In the previous Sections we introduced additional LMIs to deal with monotonicity of the PD-LF in switching and regulatory domains. Now it is possible to redefine the Feasibility Problem in Pasquini and Angel (2019) [i.e. the set of LMIs that enforce the constraints (12)-(14)] to include these modified conditions. Given the notation and discussions from the previous sections, we aim to solve the following extended feasibility problem.
Problem 1 (Extended feasibility problem) Let Σ 1 , . . . , Σ L be L extremal systems as defined in (16), and consider the associated domain partition, as described in Sect. 3. Find L Piecewise Quadratic (PWQ) functions V k (x): subject to the following constraints: Solving Problem 1 (if feasible) will give L PWQ Lyapunov functions, that when combined as in (32), return a function V λ that will be referred as a Parameter Dependent PWQ Lyapunov function for the systems in U L 1 . (42) and (44) (41), (45) and (43)]. Numerical optimization LMI solvers are used to determine a valid solution of this problem e.g. the YALMIP interface (Löfberg 2004) in MATLAB, with the SDPT3 solver.

Remark 9
The main impact in terms of computational complexity, is connected to switching domains, despite the fact that a switching domain D s enters in Problem 1 only if the associated polyhedron P D s , described in (21), is non-empty. In particular the number of switching domains increase exponentially with the size of the network and the number of thresholds, and while the operation of checking whether a polyhedron is empty or not can be performed efficiently, transform a non-empty polyhedron from its H -representation to its V -representation, as it is needed in Problem 1, is a nonpolynomial operation. We recognise this to be a potential drawback of the approach, to be addressed in further research.

Remark 10
It is reasonable to wonder if results on Lyapunov functions that are valid for hybrid systems transfer, in some way, to the case where the regulation function is smoothly close to its discontinuous counterpart (i.e. when its Hill coefficient is large enough).
Lyapunov functions with strictly negative derivative are robust to sufficiently small perturbations of the dynamics (e.g. by considering it C 0 close to PWA). Despite in our case we do not ask for strict negativity of the derivative, converse Lyapunov results, for differential inclusions (see Forni and Angeli 2017) guarantee existence of Lyapunov functions with strictly negative derivative and therefore affording some kind of robustness with respect to these perturbations. The implications of this are not explored in this context, but are certainly of interest for future research.

The set-valued derivative map
We recall and adapt the following definition from Pasquini and Angel (2019).
Definition 1 Let V k be a PWQ-LF for the extremal system Σ k as defined in (16), obtained as a solution of Problem 1. Let x ∈ D. The set-valued map if D ∈ D R , or: if D ∈ D S , withD being any regulatory domain in R(D), ∇V kD (x) defined as ∇V kD on the closure ofD, P D being the polyhedron (21) for D, F being defined as in (27) and W being the vertices set of P D 's V-representation [as expressed in (24)].

The set-valued map
• V k is connected to the derivative of V k along the system trajectories, through the following Lemma, proved in Pasquini and Angel (2019).

Lemma 1 Consider a PWQ-LF V k (x) for system (16). Letx(·) be a solution of (16) on the interval I . Then d dt V k (x(t)) is a convex combination of the elements in the set
When we consider a system σ λ ∈ C 1 L with polytopic uncertainties, Definition 1 could be naturally extended to the following: Definition 2 Let σ λ ∈ U L 1 , with U L 1 defined as in (17), let V 1 , . . . , V L be extremal LFs for the extremal systems Σ 1 , . . . , Σ L (obtained as solutions of Problem 1) and let V λ be the PD-PWQ-LF for σ λ defined in (32). Let x ∈ D. The set-valued map or: if D ∈ D S , withD being any regulatory domain in R(D), ∇V λD (x) be defined as ∇V λD on the closure ofD, P D being the polyhedron (21) for D, F being defined as in (27) and W being the vertices set of P D 's V-representation [as expressed in (24)].

Remark 11
The definition of • λ V on switching domains, takes into account the possible sliding mode directions of all the extremal systems, being P D defined as in (21). This choice is conservative, but it is legitimate as it surely contains all directions of interest, and seems necessary due to the fact that the set of sliding mode directions it is not easily described in terms of λ. Moreover we remark that for regulatory domains the set-valued map • V and the functionV are equivalent.
To complete this section we give a result about some bounds on the elements of the map • V , which will be useful in proving subsequent results.
Corollary 1 Let Σ 1 , . . ., Σ L be L extremal systems as defined in (16) and let V 1 , . . ., V L be their extremal LFs. Let the conditions of Theorem 1 be satisfied and let V λ be the . Let D ∈ D R . Then: Proof (51) follows from (40) and the fact thatV and • V are equivalent in regulatory domains.

Robust convergence properties
It is easy to extend Proposition 1 to prove that a condition like (15) holds for any system σ λ ∈ U L 1 , i.e.: However we would like to give results which take into account all the possible realizations of the uncertain system, to state convergence properties which are robust to the polytopic uncertainty.
With respect to the set-valued maps • V we introduce the following set: It is remarked that the set Ω + can be evaluated for both  (53). Consider an interval (t 0 , t 1 ), with t 0 < t 1 (t 1 could be ∞). Then: V (x(t))} for a.e. t ∈ (t 0 , t 1 ) (and more in general, for a.e. t in the whole interval of definition of the solution x(·)).

Let dV λ (x(t))
dt ≥ −ε and assume, by contradiction, that max{ For ε → 0, the set Ω + ( • λ V , ε) gives information on where the system will converge asymptotically, being in particular an outer approximation of the convergence set, and for this reason its study is of paramount importance. The same set can be defined also for the extremal Lyapunov function V k , and will be denoted with Ω + ( • k V , ε). Our goal is to characterize the set Ω + ( • λ V , ε), as λ varies in S L , in order to guarantee robust convergence result. We are interested in a qualitative analysis, therefore the following definitions are justified.
Definition 3 Consider the extremal system Σ k and its extremal Lyapunov function V k . We call D k ε the set:

Definition 4
We call D ε the union of the sets D k ε for all the extremal systems, formally: Definition 5 Consider a system σ λ ∈ U L 1 admitting a PD-PWQ-LF V λ . With analogy to Definition 3, we call D λ ε the set: Given that the dependence of Ω + ( • λ V , ε) on λ ∈ S L is extremely hard to describe, we aim to define a relation in terms of domains (i.e. the relationship between D k ε , D ε and D λ ε ). The following Theorem gives such relationship.
Theorem 2 Let Σ 1 , . . ., Σ L be L extremal systems as defined in (16)  V (x) is expressed as: in whichD is any regulatory domain adjacent to D. Being: let: Notice that ξ k, j (x) is an element of for any x ∈ D and any j ∈ {1, . . . , ν γ }.
Assume, for the sake of contradiction, that ξ k, j (x) < −ε for every k, j and x ∈ D (i.e. there is no k such that D ∈ D k ε ⊆ D ε ). Because λ ∈ S L , this yields to: which contradicts (57), proving that D ∈ D ε . We now turn to prove the same inclusion for regulatory domains.
Let D be a regulatory domain, belonging to the set D λ ε . Similarly to the previous case, it holds: implying that ∃k ∈ {1, . . . , L}, for which: This concludes the proof. Corollary 2 Given the premises of Theorem 2, it holds: The fact that: follows from the fact that any set D k ε can be obtained by selecting a vertex of the symplex S L . The fact that: directly follows from Theorem 2.
As hinted above, the set D ε contains information on the convergence set, for all the systems in U L 1 , a fact that will be formalised by Proposition 3 below, but before proceeding to state the proposition and its proof we need to introduce the concept of sink regulatory domains. In particular some regulatory domains have the property that whenever a trajectory enter them, it doesn't leave. Formally we call sink domain any regulatory domain D, The analysis inside sink domains is trivial (Casey et al. 2006) and to include them in the feasibility problem to be solved does not add any information to the study of convergence. At the same time they restrict the set of feasible extremal Lyapunov functions that satisfy the constraints in Theorem 1, as Φ(D) "moves" with the presence of uncertainties, as the ones described above, in a way that is nearly impossible to capture with a PD-PWQ-LF linear in the extremal LFs. On the other hand a bound on the position of Φ(D) could be easily computed and we know that in D the trajectories will monotonically converge to Φ(D).
For these reasons, we remove from the analysis any sink domain, including them in the convergence set, as a trajectory could always enter (or start in) a sink domain and converge to their focal point.
Proposition 3 Let Σ 1 , . . ., Σ L be L extremal systems as defined in (16) and let V 1 , . . ., V L be their extremal LFs, obtained as solution of Problem 1. Let λ ∈ S L and let V λ be the PD-PWQ-LF, as defined in (32), for σ λ ∈ U L 1 . Let D ε be the set in Definition 4 and let S D be the set of sink domains. Let x(·) be a solution of σ λ and let D(x(t)) be the function that associates the point x(t) to the domain it belongs. Then: Proof From (52) we know that: (53) for V λ . It follows easily from Lemma 1 that dV λ (x(t)) dt ∈ conv{ • λ V } for a.e. t and from Proposition 2, it follows that: where T is a set of times whose measure is finite. By definition of D ε and Theorem 2, it holds: The fact that V λ (t) is absolutely continuous (being the convex combination of absolutely continuous functions-see Pasquini and Angel (2019) implies that the derivative dV λ dt exists almost everywhere. This fact and the property that any sink domain is positively invariant, prove the thesis.

Remark 12
In the proof of Proposition 3 we used the fact that if A ⊆ [0, ∞) then: where μ is the Lebesgue measure of the set A.

Remark 13
Proposition 3 states that the trajectories will spend, asymptotically, most of the time in the domains in D ε ∪ S D, in the sense that, given a trajectory x(·), the measure of the set of times spent away from this set, shrinks to 0 as t → ∞.

Remark 14
In this work we considered polytopic uncertainties on the production rates only. However the overall analysis and LMIs framework can be easily adapted to handle uncertainties on the degradation rates only, changing the structure of the matrix δP k j D in (30), or on both the production and degradation rates, by considering all the possible combinations of their extremal values, as vertices of the uncertainty polytope.
System (59) has polytopic uncertainties on the parameters b 2 and b 3 , so there are four extremal systems, corresponding to the possible combinations of b 2 and b 3 extrema.
Problem 1 is set up and hence the following constraints will be asked:  (42) and (43). In this example, these domains are: 3. Continuity of the extremal Lyapunov functions along the cycles in the STG of the network, which converts into asking continuity on the walls connecting the domains D 5 , D 6 , D 7 and D 8 , relatively to the partition in Fig. 7, as per constraints (44).
We can notice that the regulatory domain: is a sink domain, and for this reason it is removed from the formulation of Problem 1. By solving Problem 1, it is possible to find four extremal Lyapunov functions, from which we can construct a PD-LF V λ . In Fig. 8  realizations of the uncertainty and in Fig. 9 it is shown how the parameter dependent Lyapunov function V λ , is eventually non-increasing along these trajectories. More numerical details on the LMI formulation and potential solution are given in the "Appendix" Choosing ε = 10 −7 , it is possible to find the sets D k ε , for all the extremal Lyapunov functions V k , and consequently we can determine the set D ε , corresponding to: where: Given the cyclic nature of the system trajectories around D s 2 , and because • k V (x), for each extremal LF, approach 0 as x tends to D s 2 , the set D ε contains also the regulatory domains adjacent to D s 2 (i.e. {D 5 , D 6 , D 7 , D 8 }).
However from the proof of Theorem 2, it follows: so that, by studying the set where: we can infer that the solutions always converge, in the sense of measure, to the domain D s 2 . This is also a consequence of the fact that • k V is described, in the regulatory domains around D s 2 , by non-singular quadratic forms, and it is empty in the boundaries connecting these regulatory domains. To give a visual idea of the result, in Fig. 10, different level sets of the function max with decreasing values of ε. The projection on (x 1 , x 2 ) is considered, as V λ does not depend on x 3 in those particular regulatory domains. This analysis is supported by the trajectories represented in Fig. 11, where the highlighted cyan domain is D s 2 , to which we can see that most of the trajectories converge. It should be noticed that, for the same initial condition x 0 = 3 1.1 0.5 T , the system may converge to the switching domain D s 2 or to the sink domain D 2 , depending on the uncertainty realization.

Remark 15
In Casey et al. (2006), a conjecture-proved true in Wang and Wang (2013) -is given with sufficient conditions for the weak asymptotic stability of a switching domain D s , under the assumption that no cycles are present in the STG reduced to the neighbours of such switching domain. In the case of the above example however, the domain D s 2 has two cycles involving neighbour domains, hence the result of Casey et al. (2006) and Wang and Wang (2013) cannot be applied.

Conclusion
Due to context-dependence and measurement limitations, models of biological systems are affected by parameter uncertainties and different kind of disturbances. In this chapter we considered a PWA model of genetic regulatory networks, and assumed the presence of polytopic uncertainties affecting the production rates. The LMI framework of Pasquini and Angel (2019) has been extended, allowing to determine a Parameter Dependent Piecewise Quadratic Lyapunov function, for the whole polytope of uncertain parameters. With this function we can conclude on the convergence properties of the trajectories to a particular set of domains, independently from the values assumed by the perturbed parameters. This can rule out a big chunk of the positive orthant from the convergence set, with further analysis that can be done in the resulting convergence set, to better characterise the system dynamics. In fact, as has been done in the double feedback example, the study of the level sets for the derivative map can allow to conclude more precise results on the system convergence. These results can guide both in the synthesis of de-novo genetic circuits and in predicting the behavior of natural ones.
However there are a couple of issues in the analysis, that deserve to be addressed in future research. Firstly, as already addressed in Remark 3, Assumption 1 considers that the thresholds are unchanged among the extremal systems, which we recognize to be unrealistic. An extension of the framework, to allow polytopic uncertainties in the thresholds-e.g. by giving a lower and an upper-bound on each θ -would certainly be of interest. However this would require a significant structural change of the framework, as continuity constraints would need to be applied on an undefined threshold, with their convexification leading to non-linear conditions, moving away from an LMI setup. If one chooses to drop any continuity constraint instead, the entire theoretical scaffolding should be changed, as infinite jumps of the Lyapunov functions would now be allowed, while at the same time, entry and exit points of sliding modes on the surface of discontinuity should be taken into account. These matters should be investigated in future research.
Secondly the constraints relative to the monotonicity of the Lyapunov function along sliding mode can be very conservative, and extremely expensive even to define. This is due to the fact that it is difficult to characterise the set of sliding directions in terms of the uncertain parameters, hence we need to impose this condition on all possible sliding directions, for all extremal systems. A better characterisation of this set-and the consequent relaxation of the LMI conditions-would increase the set of feasible solutions, while decreasing the computational cost of the approach.
In connection with this latter point, and Remark 9, general efficiency of the framework implementation should be pursued as many Genetic Regulatory Networks in nature are composed of tens of nodes. The framework described here can easily handle networks of up to six or seven nodes, but the computational cost becomes prohibitive after that. This work should be in fact considered as a proof of concept, and we believe that there are many possibilities to increase its efficiency.
Ultimately, future works should also consider control applications of the analytical results given above. Of particular interest is the problem of in-vivo control of Genetic Regulatory Networks, namely the addition, to the original GRN, of a second network whose goal is to obtain a desired overall behavior (e.g. control of an unstable steady state). With the current framework one could attempt a trial and error approach, meaning that a structure for the control-network can be chosen, the overall convergence properties studied through the described Lyapunov approach, and the control network modified accordingly, until the desired behavior is reached. We recognize that this approach is extremely inefficient and the possibility of embedding the control problem into a complementary set of LMIs should be addressed in future research.
Funding Open access funding provided by Royal Institute of Technology. This work was supported in part by the Engineering and Physical Sciences Research Council (EPSRC), under award 1895642, and in part by the Department of Electrical and Electronic Engineering at Imperial college of London.

Conflicts of interest The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
where Γ D 1 is the ray matrix of the homogenization cone C D 1 , which can be constructed as in Eq. (3), in this case: for which P 1 D 1 and d 1 D 1 are decision variables, while: is another decision variable of the problem, with the only constraints that its elements are all non-negative. An LMI interpreter (such as Yalmip Löfberg 2004 in MATLAB) with a semidefinite programming solver (such as SDPT3 Toh et al. 1999) can be used to impose this kind of constraints and search for a solution. In regards to (45), consider for example the crossed conditions between the first and the second extremal Lyapunov functions (i.e. V 1 and V 2 ), for the regulatory domain D 4 . This converts to: where: where, as previously, P 2 D 4 , P 1 D 4 , d 1 D 4 , d 2 D 4 and M 12 D 4 are decision variables, while: We remark that conditions (68) and (69) are asked for all regulatory domains (except D 2 ) and all extremal Lyapunov functions (with the crossed conditions, asked for any possible couple of extremal Lyapunov functions). For switching domains the first step is to check if, for any given D S , the polyhedron P D S described in Eq. (21) is non-empty. This operation can be done efficiently using available software libraries for polyhedra (e.g. MPT3 Herceg et al. 2013). In this case the only two switching domains with such properties are: The H -representation (21) of the polyhedron P D S , is then converted to its V -representation, and again this operation can be performed using dedicated software toolboxes (e.g. MPT3). Condition (43) is then asked, for all vertices in such V -representation and all extremal systems. It is remarked that the polyhedron P D S is independent of the particular extremal Lyapunov function, by construction.
Moreover for the two switching domains, continuity of the extremal Lyapunov functions is asked. Consider for example the switching domain D S 2 and the first extremal LF. Condition (42) converts to: where: while P 1 D is the matrix:  (70) is then asked for all extremal Lyapunov functions and for all the other switching domains with non-empty P D S , in this case only D S 1 . Ultimately the conditions on continuity along cycles of the STG should be asked. Given the STG (which can be obtained through the algorithm explained in De Jong et al. 2004), we can search for all the cycles in it using efficient algorithms (see for example Tiernan 1970), and for each cycle ask for the extremal Lyapunov functions to be continuous on the switching domains which will eventually be crossed. In this case this converts to asking for the extremal Lyapunov function to be continuous on the walls connecting the regulatory domains D 5 , D 6 , D 7 and D 8 . Consider for example the wall between D 5 and D 6 , for the extremal Lyapunov function V 1 . The continuity condition becomes: where D 56 is the wall between D 5 and D 6 , while: This type of condition is repeated for the remaining walls and for all other extremal Lyapunov functions. Once all the LMIs and matrix equalities are collected, the feasibility problem can be solved (with the aforementioned software) and a description of the parameter dependent Lyapunov function (in terms of the variables P k D , d k D and ! k D for all regulatory domains D, for all extremal Lyapunov functions V k ). As an example, below the description of the first extremal Lyapunov function: