On the robust existence of piecewise quadratic Lyapunov functions for hybrid models of gene regulatory networks

TL;DR

This paper develops conditions for the existence of piecewise quadratic Lyapunov functions to analyze the stability of uncertain hybrid gene regulatory network models, accounting for parameter uncertainties and sliding modes.

Contribution

It provides a novel theoretical framework for stability analysis of uncertain hybrid gene regulatory networks using piecewise quadratic Lyapunov functions.

Findings

Conditions for Lyapunov function existence are derived for all parameter uncertainties.

The approach accounts for sliding modes in hybrid systems.

An example demonstrates the applicability of the theoretical results.

Abstract

In this work we addressed the problem of stability analysis for an uncertain piecewise affine model of a genetic regulatory network. In particular we considered polytopic parameter uncertainties on the proteins production rate functions, giving conditions for the existence of a piecewise quadratic Lyapunov function for any possible realization of the system. In the spirit of other works in literature, the resulting conditions will be given on the vertices of the parameter polytope, while still taking into consideration the piecewise nature of the Lyapunov function and the presence, in general, of sliding modes solutions. An example is shown to prove the validity and applicability of the theoretical results.

Equations

$${\mathbb{R}}_{+}^n:=\{v\in {\mathbb{R}}^n|v\ge 0\}$$

$$cone\{\Omega \}:=\{v\in {\mathbb{R}}^n|\exists \alpha \in {\mathbb{R}}_{+}^n\text{ s.t. }v=\sum _{i=1}^n{\alpha }_i{\omega }_i\}$$

$$S_r:=\{v\in {\mathbb{R}}_{+}^r|\sum _{i=1}^r{v}_i=1\}$$

$$conv\{\Omega \}:=\{v\in {\mathbb{R}}^n|\exists \alpha \in S_n\text{ s.t. }v=\sum _{i=1}^n{\alpha }_i{\omega }_i\}$$

$$x^{T}Mx>0(x^{T}Mx\ge 0)\forall x\in {\mathbb{R}}^n$$

$$P:=\{x\in {\mathbb{R}}^n|Ax\le b,A\in {\mathbb{R}}^{m\times n},b\in {\mathbb{R}}^m\}$$

$$P:=conv\{V\}+cone\{R\}$$

$${\dot{x}}_i=f_i(x)-c_i{x}_i,x_i\in {\mathbb{R}}_{+}$$

$$f_i(x):=\sum _{j=1}^{l_i}{b}_{i,j}\cdot s_{i,j}(x)$$

$$s^{+}(x_{i},\theta_{i,k})=\begin{cases}\begin{array}[]{lll}1,&\text{ if }&x_{i}>\theta_{i,k}\\ 0,&\text{ if }&x_{i}<\theta_{i,k}\end{array}\end{cases}$$

$$s^{-}(x_i,{\theta }_{i,k})=1-s^{+}(x_i,{\theta }_{i,k})$$

$$\dot{x}=f(x)-Cx,x\in {\mathbb{R}}_{+}^n$$

$$f(x)=\left[\begin{array}{c}{f}_1(x)\\ ⋮\\ f_n(x)\end{array}\right]$$

$$\dot{x}=f_D-Cx,x\in {\mathbb{R}}_{+}^n$$

$$\dot{x}\in H(x)$$

$$H(x)=\begin{cases}\begin{array}[]{lll}\{f_{D}-Cx\},&\text{if }D\in\mathcal{D}_{R}\\ conv\{f_{D^{\prime}}-Cx\ |\ D^{\prime}\in R(D)\},&\text{if }D\in\mathcal{D}_{S}\end{array}\end{cases}$$

$${\Sigma }^k:\dot{x}=f^k(x)-Cx,x\in {\mathbb{R}}_{+}^n$$

$${\mathcal{C}}_1^L:=\{{\sigma }^{\lambda }:\dot{x}=f^{\lambda }-Cx,x\in {\mathbb{R}}_{+}^n\}$$

$$f^{\lambda }:=\sum _{k=1}^L{\lambda }_k{f}^k(x),\lambda \in S_L$$

$$\dot{x}=\sum _{k=1}^L{\lambda }_k{f}_D^k-Cx$$

$$f=\sum _{D^{\prime }\in R(D_s)}{\alpha }_{D^{\prime }}\left(\sum _{k=1}^L{\lambda }_k{f}_{D^{\prime }}^k\right)-Cx$$

$$f=\sum _{D^{\prime }\in R(D_s)}\sum _{k=1}^L{\alpha }_{D^{\prime }}{\lambda }_k{f}_{D^{\prime }}^k-Cx$$

$$\begin{array}{r}\sum _{D^{\prime }\in R(D_s)}\sum _{k=1}^L{\alpha }_{D^{\prime }}{\lambda }_k=\sum _{D^{\prime }\in R(D_s)}{\alpha }_{D^{\prime }}\sum _{k=1}^L{\lambda }_k=\sum _{D^{\prime }\in R(D_s)}{\alpha }_{D^{\prime }}=1\end{array}$$

$$f=\sum _{D^{\prime }\in R(D_s)}\sum _{k=1}^L{\gamma }_{D^{\prime }}^k{f}_{D^{\prime }}^k-Cx$$

$$P_{\gamma,D_{s}}:=\left\{\gamma\in\mathbb{R}^{L\cdot q}:\begin{cases}\begin{array}[]{lcr}\begin{bmatrix}\overline{F^{1}}&\overline{F^{2}}&\dots&\overline{F^{L}}\end{bmatrix}\gamma&=&\begin{bmatrix}\bar{c}\end{bmatrix}\\ 1^{T}\gamma&=&1\\ \gamma&\geq&0\\ \end{array}\!\!\!\!\!\!\!\!\end{cases}\right\}$$

$$[\bar{c}{]}_i=c_i{\theta }_{i,k},i\in I_{D_s}$$

$$\overline{F^k}=\left[\begin{array}{ccc}(f_{D_{i_1}}^k{)}^{I_{D_s}}& \dots & (f_{D_{i_q}}^k{)}^{I_s}\end{array}\right]$$

$$P_{\gamma ,D_s}:=conv\{w_1,w_2,\dots ,w_{{\nu }_{\gamma }}\}$$

$$H^{\lambda }(x)\subseteq conv\{H^1(x),H^2(x),\dots ,H^L(x)\}$$

$$\varphi (D):=C^{-1}{f}_D$$

Full text

Copyright ©2019 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to servers or lists, or reuse of any copyrighted component of this work in other works.

I INTRODUCTION

Due to surrounding environment and measurement noises, biological models often contain parameters that are uncertain and subject to variability [1, 2], and for this reason it is appropriate to investigate the stability properties of such systems in a robust sense.

One of the most important concept in systems biology is that of Genetic Regulatory Networks (GRNs), as these networks describe the regulatory interactions between genes and proteins in a cellular environment [3].

In [4] and [5] the authors, using a piecewise affine (PWA) model of the dynamics of a GRN, employed a state transition graph (STG) – which is unchanged in a large set of parameters – to qualitatively describe the system trajectories. They found conditions that allow to determine stability properties of equilibria, which are related only to the STG.

In [6], we considered this PWA model and, using some interesting mathematical tools presented in [7, 8], we built an LMI framework that, solving a feasibility problem, gave the description of a piecewise quadratic (PWQ) Lyapunov function for the system. Nevertheless the solution found was strictly dependent on the system parameters, as these determine the values of the matrices defining the LMIs.

Many authors found ways to extend the Lyapunov approach, in order to deal with robustness stability when the system is subject to polytopic uncertainties. In [9] the authors considered a linear system, in which the system matrix belongs to a polytope of matrices. They gave a set of LMI conditions that, if verified, return the description of a Lyapunov function which is linearly dependent from the uncertain parameters, proving the Hurwitz stability of the entire set of system matrices in the aforementioned polytope. In [10] and [11] the authors, considering the same kind of polytopic uncertainties, extended the approach considering Lyapunov functions that are dependent in an homogeneous polynomial way on the parameters. Despite considering linear systems with uncertain system matrices, these works did not consider switched systems and the problems associated with those, such as: the presence of sliding modes solutions and the continuity constraint to be asked in certain regions for a (generally) discontinuous piecewise quadratic Lyapunov function. In the following we will consider polytopic uncertainties on the production rate functions, in a PWA model of an N proteins system, and we present two approaches, in the spirit of [9, 10, 11], to prove the robust existence of a PWQ-LF under the aforementioned uncertainties, extending the LMI framework we described in [6].

The work is organized as follows: in Section II some mathematical background and notation are presented, while in Section III the piecewise affine model of the dynamics of a genetic regulatory network is briefly recalled. In Section IV the problem is formulated in a mathematically precise way and in Sections V and VI two solutions to the problem are illustrated. Finally Section VII shows the application of the results to a numerical example while Section VIII concludes the paper highlighting some possible further developments.

II NOTATION AND MATHEMATICAL BACKGROUND

Let $v=\begin{bmatrix}v_{1}&\dots&v_{n}\end{bmatrix}^{T}$ be a vector in $\mathbb{R}^{n}$. The notation $v\geq 0$ means that $v_{i}\geq 0,\ \ \forall i\in\{1,\dots,n\}$. With $\mathbb{R}^{n}_{+}\subseteq\mathbb{R}^{n}$ we denote the set:

[TABLE]

Let $\Omega:=\{\omega_{1},\dots,\omega_{n}\}\subset\mathbb{R}^{n}$ be a finite set of vectors in $\mathbb{R}^{n}$. The conic hull of $\Omega$ is defined as:

[TABLE]

With $S_{r}$ we denote the standard simplex of dimension $r$, namely the set:

[TABLE]

The convex hull of $\Omega$ is defined as:

[TABLE]

For any set $D$: $int(D)$ indicates the interior of $D$, $\partial D$ indicates its boundary and $cl(D)=int(D)\cup\partial D$ its closure. The notations $cl(D)$ and $\overline{D}$ can be used interchangeably.

Let $M\in\mathbb{R}^{n\times n}$ be a symmetric matrix of size $n$. $M$ is said to be positive definite (semidefinite), and indicated with $M\succ 0$ ($M\succeq 0$), if:

[TABLE]

$M$ is said to be negative definite (semidefinite), and indicated with $M\prec 0$ ($M\preceq 0$), if $-M$ is positive definite (semidefinite).

A polyhedron $P$ in $\mathbb{R}^{n}$ is the set of solutions of a system of linear inequalities, formally:

[TABLE]

The description (6) is called $H$-representation of $P$ and can always be converted in a $V$-representation [12, 13], formally:

[TABLE]

where $V$ is the set of vertices of $P$ and $R$ the set of its rays. A procedure exist [7] to associate an higher dimension cone, called homogenization cone, to any polyhedron $P$ and such cone will be used in the following to give conewise condition on quadratic functions.

III PIECEWISE AFFINE MODEL OF A GRN

*In this section we briefly recall the piecewise affine model considered for a genetic regulatory network. Most of this section is taken from [4, 6], to which the reader is referred for a deeper discussion.

*Consider a system of $n$ proteins and, for the single protein $X_{i}$, consider its concentration $x_{i}$ to evolve according to the dynamics:

[TABLE]

in which $x=\begin{bmatrix}x_{1}&x_{2}&\dots&x_{n}\end{bmatrix}^{T}$ is the vector of proteins concentrations, $f_{i}(x)\geq 0$ describes the production rate of $X_{i}$ and $c_{i}>0$ is its degradation rate, that takes in consideration both diluition inside the cell and natural degradation of $X_{i}$. We assume $f_{i}(x)$ to be of the form:

[TABLE]

for which $s_{i,j}(x)$ is a sum of products of step functions:

[TABLE]

[TABLE]

When considering all the $n$ proteins, we have the following model:

[TABLE]

in which $C\in\mathbb{R}^{n\times n}$ is a diagonal matrix with strictly positive diagonal entries and $f(x):\mathbb{R}^{n}_{+}\to\mathbb{R}^{n}_{+}$ is defined as:

[TABLE]

The structure of (9) and the step functions (10) and (11), naturally partitions the state space in boxes, defined by the thresholds $\theta$. Such boxes (possibly with empty interior) will be called domains. If we consider a domain $D$, in which none of the $x_{i}$s assume a threshold value, such domain will be called regulatory, while if some $x_{i}$s assume a threshold value, the domain $D$ will be called switching. The set of regulatory domains is indicated with $\mathcal{D}_{R}$, while the set of switching domains is indicated with $\mathcal{D}_{S}$. For a domain $D_{s}\in\mathcal{D}_{S}$, we denote with $I_{D_{s}}$ the set of indexes of the $x_{i}$s that assume a threshold value in $D_{s}$.

Inside a particular regulatory domain $D$ the function $f(x)$ is a constant vector and will be indicated with $f_{D}$, giving rise to the regulatory dynamics:

[TABLE]

if $x\in D$. However in a switching domain the function $f(x)$ is not properly defined in at least one of its components, so the system (12) needs to be extended to a differential inclusion [14], namely:

[TABLE]

with:

[TABLE]

where $R(D)$ indicates the set of regulatory domains adjacent to the switching domain $D$. In [5, 15, 4] the possibility to associate a State Transition Graph (STG) to the system is described. An STG is a graph in which: the nodes are associated to the domains of (12), while every edge is associated to the existence of a solution in the sense of Filippov [14], which connects two domains, without crossing any other domain in between.

In [6] we developed an LMI framework that consents to find a Piecewise Quadratic (PWQ) Lyapunov function $V(x)$ for system (12). Such method is based on enforcing a set of LMIs and Matrix Equalities (MEs) constraints regarding: the continuity (only where required, in general the function will be discontinuous) and the monotonicity of $V$.

However, the feasibility of the problem is strictly related to the system parameters that, as stated at the beginning of this paper, are mostly uncertain and possibly varying with time. The next sections will deal with assuring the feasibility of the LMIs framework when the system is subject to polytopic uncertainties of the production rate function.

IV PROBLEM FORMULATION

Let $C\in\mathbb{R}^{n\times n}$ be a diagonal matrix with strictly positive diagonal entries. Let $f^{1}(x)$, $\dots$, $f^{L}(x)$ be $L$ different production rate functions of the type (13). Let $\Sigma^{k}$ be the system:

[TABLE]

We call $\Sigma^{k}$ an extremal system. Let $\mathcal{C}_{1}^{L}$ be the set of systems:

[TABLE]

for which $f^{\lambda}$ has the form:

[TABLE]

in which $S_{L}$ is the standard simplex of dimension $L$. Given the domain dependent nature of $f^{k}(x)$, any system $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$ has the dynamics:

[TABLE]

inside the regulatory domain $D$. We make the following assumption on $\mathcal{C}_{1}^{L}$.

Assumption 1**.**

All the systems in $\mathcal{C}_{1}^{L}$ have the same state transition graph, the same thresholds and the same domains.$\hfill\blacksquare$

Our goal is to solve the following:

Problem 1**.**

Find conditions under which any system in $\mathcal{C}_{1}^{L}$ admits a Piecewise Quadratic Lyapunov function. $\hfill\blacksquare$

First we need to discuss how to describe sliding modes for a generic $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$, then we will explain two approaches to solve Problem 1, giving also a description of the resulting Lyapunov function.

IV-A Sliding modes for systems in $\mathcal{C}_{1}^{L}$

Let $\sigma^{\lambda}$ be a system in $\mathcal{C}_{1}^{L}$. We know that if $x\in D_{S}$, with $D_{S}$ being a switching domain, then $\sigma^{\lambda}$ should be extended to a differential inclusion [14], and $\dot{x}\in H^{\lambda}(x)$, where $H^{\lambda}(x)$ is the set of directions $f$ satisfying the following equation:

[TABLE]

with $\alpha\in S_{|R(D_{s})|}$ and $\lambda\in S_{L}$. Equation (21) can be written as:

[TABLE]

and, given that both $\alpha$ and $\lambda$ belong to standard simplices, we have:

[TABLE]

so that (21) can be expressed as:

[TABLE]

with $\gamma\in S_{L\cdot q}$ and $q:=|R(D_{s})|$. Following an approach similar to [6], the sliding mode directions of $\sigma^{\lambda}$ are certainly included in the ones described by (24), in which $\gamma$ belongs to the following polyhedron:

[TABLE]

where $\bar{c}$ is a vector, the i-th component of which is:

[TABLE]

where $c_{i}$ is the degradation rate of the protein $X_{i}$ and $\theta_{i,k}$ is the threshold value assumed by $x_{i}$ in the switching domain $D_{s}$, while $\overline{F^{k}}$ is the matrix:

[TABLE]

in which $(f^{k}_{D})^{I_{D_{s}}}$ is the vector obtained by selecting only the components from $(f^{k}_{D})$, indexed by the set of switching variables $I_{D_{s}}$ and $\{D_{i_{1}},\dots,D_{i_{q}}\}$ is the set of regulatory domains adjacent to $D_{S}$. Because any $H$-representation of a polyhedron $P$ can be converted to its $V$-representation [12, 7],and given the boundedness of $P_{\gamma,D_{s}}$, we get:

[TABLE]

Remark 1**.**

A direct consequence of the above analysis is that:

[TABLE]

meaning that any element in $H^{\lambda}(x)$ is a convex combination of elements in the sets $H^{1}(x),H^{2}(x),\dots,H^{L}(x)$, where $H^{k}(x)$ is the set (16) for the extremal system $\Sigma^{k}$. Nevertheless this does not guarantee that any sliding mode direction of $\sigma^{\lambda}$ is a convex combination of the sliding mode directions of the extremal systems $\Sigma^{1}$, $\dots$, $\Sigma^{L}$ (i.e. not all the directions in $H^{k}(x)$ are sliding mode directions for $\Sigma^{k}$). $\hfill\blacksquare$

V Common Lyapunov function approach

Intuitively if a common piecewise quadratic Lyapunov function $V$ exists for all the extremal systems $\Sigma^{k},\ k\in\{1,\dots,L\}$, then $V$ is a possible Lyapunov function for any $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$. Unluckily, just searching for a common Lyapunov function will, almost certainly, result in an unfeasible problem, due to the fact that different systems in $\mathcal{C}_{1}^{L}$ will have different equilibria. The following definition is instrumental to bypass this problem:

Definition 1**.**

A regulatory domain $D$ is said to be a sink domain if $\phi(D)\in D$.$\hfill\blacksquare$

We recall from [5, 4] that the focal point $\phi(D)$ of the regulatory domain $D$ for (14) is defined as:

[TABLE]

Sink domains are the ones that drastically reduce the possibility of finding a common Lyapunov function for the extremal systems. On the other hand the dynamics inside sink domains is well characterized [4], with the property that any trajectory entering a sink domain will not leave it. This means that we don’t need to actually define a Lyapunov function here, as we already know how the system behaves. It is a direct consequence of Assumption 1 that the set of sink domains is the same for all the systems in $\mathcal{C}_{1}^{L}$. The approach is then to search a piecewise quadratic function $V$, being a common Lyapunov function for the extremal systems, restricted only to non-sink domains.

Suppose that such function $V$ exists. We need to prove that $V$ is in fact a Lyapunov function for any $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$. Inside a regulatory domain $D$ we have:

[TABLE]

With respect to switching domains, by Remark 1, it is not sufficient for $V$ to be non-increasing along any sliding mode direction of any of the extremal systems, as some sliding mode direction of $\sigma^{\lambda}$ could be not included in those. Let $D_{S}$ be a switching domain, for which $P_{\gamma,D_{S}}$ is non empty, and let $L_{D_{s},j}$ be the matrix:

[TABLE]

in which $F$ is:

[TABLE]

and $F^{k}$ is the matrix:

[TABLE]

for the extremal system $\Sigma^{k}$. If, together with the continuity of $V$ in $D_{S}$, the following set of LMIs is satisfied [7, 6]:

[TABLE]

in which $\Gamma_{D_{s}}$ is the ray matrix of the homogenization cone of $D_{s}$ and $M_{D_{s},j}$ is any entrywise non-negative and symmetric matrix, then $V$ is non-increasing along any direction corresponding to $\gamma\in P_{\gamma,D_{s}}$ and because $P_{\gamma,D_{s}}$ contains the $\gamma$s that give rise to all the sliding mode directions of $\sigma^{\lambda}$ in $D_{s}$, this implies that $V$ is non-increasing along any sliding mode trajectories of $\sigma^{\lambda}$ in $D_{S}$.

Lastly continuity constraints on $V$ are dictated by the STG of the system, which is unchanged by Assumption 1, and by the switching domains $D_{S}$ for which the polyhedron $P_{\gamma,D_{S}}$ is non empty. Trivially, because of Assumption 1 and the above discussion on sliding modes, such $V$ will respect all the required continuity constraints. The so found $V$ will then be a common piecewise quadratic Lyapunov function between all the systems in $\mathcal{C}_{1}^{L}$.

VI Additional LMIs approach

Despite the simplicity of the approach in Section V, the constraint of $V$ being in common between the extremal systems, can be too restrictive. Hence a second approach consists in finding $L$ different piecewise quadratic Lyapunov functions $V^{1}$, $V^{2}$, $\dots$, $V^{L}$ for the extremal systems that, together with the usual constraints described in [6], will satisfy a set of additional constraints, guaranteeing the existence of a function $V^{\lambda}$ for any $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$.

This approach is a generalization of the first one in the sense that if a solution exists to the problem defined in Section V, then it will also be a feasible solution of the problem that we are now defining.

The function $V^{\lambda}$ is imposed to be a conic combination of $V^{1}$, $V^{2}$, $\dots$, $V^{L}$, formally:

[TABLE]

We call the functions $V^{1}$, $V^{2}$, $\dots$, $V^{L}$ extremal Lyapunov functions relative to the extremal systems $\Sigma^{1},\Sigma^{2},\dots,\Sigma^{L}$. For clarity of the following explanation, the extremal Lyapunov function relative to the extremal system $\Sigma^{k}$ will be described by:

[TABLE]

In the subsequent discussion we are going to express all the constraints on $V^{\lambda}$ in terms of the extremal Lyapunov functions.

VI-A Non-increasing in switching domains

The first set of constraints is relative to sliding mode solutions. In particular we need to ensure that the Lyapunov function $V^{\lambda}$ of $\sigma^{\lambda}$, is non-increasing along any sliding mode direction. To guarantee this we use a technique similar to the one applied in Section V, asking that, in the switching domain $D_{S}$, all the extremal Lyapunov functions $V^{k}$ are non-increasing along any direction:

[TABLE]

in which F is the one described in (33) and $\gamma\in P_{\gamma,D_{S}}$, with $P_{\gamma,D_{s}}$ described by (25). If the following set of LMIs is satisfied:

[TABLE]

with:

[TABLE]

then the desired condition is satisfied as well.

VI-B Continuity constraints

As discussed in Section V, being the STG unchanged, and because of the constraints of Section VI-A, the required continuity constraints are the same for all the systems in $\mathcal{C}_{1}^{L}$. For such reason all the functions $V^{k}$ are continuous in the same domains, and being the function $V^{\lambda}$ a conic combination of the functions $V^{k}$, it will be continuous in the same domains as well, proving that however $\eta_{1}$, $\eta_{2}$, $\dots$, $\eta_{L}$ are chosen in (36), all the required continuity constraints hold for $V^{\lambda}$.

VI-C Non-increasing in regulatory domains

The last set of constraints is related to the monotonicity of $V^{\lambda}$ inside regulatory domains. Let $D$ be a regulatory domain for the systems in $\mathcal{C}_{1}^{L}$. The extremal Lyapunov functions satisfy [6]:

[TABLE]

where $\tilde{P}_{D}^{k}$ is the matrix:

[TABLE]

and $M_{D}^{k}$ is a symmetric and entrywise-non-negative matrix. Our goal is to give a set of constraints which guarantees that, for any $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$, it exists $\eta\in\mathbb{R}^{L}_{+}$ such that, for $V^{\lambda}$ as expressed in (36), the following LMI is satisfied:

[TABLE]

for any regulatory domain $D$, in which:

[TABLE]

and:

[TABLE]

In order to do this the matrix $\tilde{P}_{D}^{\lambda}$ should be expressed in terms of the matrices $\tilde{P}_{D}^{k}$s.

Let $\left[\tilde{P}_{D}^{\lambda}\right]_{ij}$ denotes the block in position $(i,j)$ of $\tilde{P}_{D}^{\lambda}$ so that:

[TABLE]

If, respectively, $\left[\tilde{P}_{D}^{k}\right]_{ij}$ denotes the block in position $(i,j)$ of $\tilde{P}_{D}^{k}$, it is immediately clear that:

[TABLE]

For the other terms more computations are needed. In particular for $\left[\tilde{P}_{D}^{\lambda}\right]_{12}$ it holds:

[TABLE]

Knowing that $\lambda\in S_{L}$ we can rewrite (48) as:

[TABLE]

Rearranging the terms of (49) we get:

[TABLE]

Defining:

[TABLE]

we obtain:

[TABLE]

Considering the following properties of $\delta f^{kj}$:

[TABLE]

after a few manipulations of (52) we obtain:

[TABLE]

Now assume that $\lambda\in int(S_{L})$ (i.e. $\lambda_{k}\neq 0,\forall k\in\{1,\dots,L\}$). If we choose $\eta_{k}$ as:

[TABLE]

substituting (55) in (54), and defining:

[TABLE]

gives:

[TABLE]

in which:

[TABLE]

Moreover defining $\delta d_{D}^{kj}=d_{D}^{j}-d_{D}^{k}$, with similar reasoning we get:

[TABLE]

If we put together (47), (57) and (59), we can rewrite $P_{D}^{\lambda}$ as:

[TABLE]

in which:

[TABLE]

Using (60), condition (43) can be written as:

[TABLE]

and rearranging the terms in (62) we obtain:

[TABLE]

Because constraints (41) hold and because any $\xi_{kj}>0$, then satisfying the following set of LMIs:

[TABLE]

with $M^{kj}$ being an entrywise non-negative and symmetric matrix, is sufficient to guarantee that (43) is satisfied.

Remark 2**.**

Being:

[TABLE]

the Lyapunov function $V^{\lambda}$ for $\sigma^{\lambda}$ is described by:

[TABLE]

Considering $\lambda\in int(S_{L})$, and considering that $\alpha V^{\lambda}$ is still a Lyapunov function for $\sigma^{\lambda}$, $\forall\alpha>0$, we have that:

[TABLE]

is a Lyapunov function for $\sigma^{\lambda}$. $\hfill\blacksquare$

Remark 3**.**

The analysis done in the Section VI has been done for $\lambda\in int(S_{L})$. Anyway if $\lambda\in\partial S_{L}$ it is possible to prove that the problem reduces to a lower dimensional one – meaning that if, for example, one of the $\lambda_{k}=0$, the problem is equivalent to a problem with $L-1$ extremal systems instead of L – and the relative constraints are contained in the ones already asked. $\hfill\blacksquare$

Remark 4**.**

As in Section V, if the constraints relative to sink domains are dropped, the possibility to find a feasible solution drastically improves. $\hfill\blacksquare$

VII EXAMPLE

VII-A System with sliding mode

Consider the system:

[TABLE]

where:

[TABLE]

We want to guarantee the existence of a Lyapunov function for $\lambda\in S_{4}$. We apply the approach of Section VI so, given a system $\sigma^{\lambda}\in\mathcal{C}_{1}^{L}$, we define a set of additional constraints on the extremal Lyapunov functions, guaranteeing that their convex combination with weights $\lambda_{k}$ is a Lyapunov function for $\sigma^{\lambda}$. The algorithm is able to find a feasible solution in the extremal Lyapunov functions showed in Figures 1–4.

For the system $\sigma^{\lambda}$, the function:

[TABLE]

is then a Lyapunov function. To test the validity of the results we generate values of $\lambda\in S_{4}$ and evaluate the function $V^{\lambda}$ (Figure 6) along the trajectories of system (68) (Figure 5).

VIII CONCLUSIONS AND FUTURE WORK

We considered a piecewise affine model of the dynamics of a GRN, subject to polytopic uncertainties of the production rate function. We addressed the problem of finding conditions that guarantee the existence of a piecewise quadratic Lyapunov function for any of the systems in the set $\mathcal{C}_{1}^{L}$, the dynamics of which is obtained through a convex combination of the dynamics of a finite number of, what we called, extremal systems. We proposed two solutions to this problem: the first consists in searching for a common Lyapunov function among the extremal systems, limited only to particular domains, while the second approach is a generalization of the first one and consists in searching for different Lyapunov functions for all the extremal systems that, when linearly combined, can generate Lyapunov functions for all the systems in $\mathcal{C}_{1}^{L}$. To conclude we proposed a numerical example, to show the validity of the described approach.

Future works should focus on the properties shared by the Lyapunov functions for the systems in $\mathcal{C}_{1}^{L}$, aiming to prove a result of robust convergence for the system.