B Bedessem

A Stéphanou

A mathematical model of HIF-1α-mediated response to hypoxia on the G1/S transition

Keywords: Hypoxia, HiF-1α, G1/S transition, quiescence, cell cycle arrest, cancer

.

The HiF-1α factor is the key actor of the intracellular hypoxia pathway.

As tumor cells undergo chronic hypoxic conditions, HiF-1α is present in higher level in cancer than in normal cells. Besides, it was shown that genetic mutations promoting overstabilization of HiF-1α are a feature of various types of cancers [START_REF] Brdos | Negative and positive regulation of hif-1: a complex network[END_REF]. Finally, it is suggested that the intracellular level of HiF-1α can be related to the aggressiveness of the tumors [START_REF] Zhong | Overexpression of hypoxia-inducible factor 1alpha in common human cancers and their metastases[END_REF][START_REF] Hashimoto | Hypoxia induces tumor aggressiveness and the expansion of cd133-positive cells in a hypoxia-inducible factor-1-dependent manner in pancreatic cancer cells[END_REF][START_REF] Bao | Hypoxia-induced aggressiveness of pancreatic cancer cells is due to increased expression of vegf, il-6 and mir-21, which can be attenuated by cdf treatment[END_REF][START_REF] Chen | Hypoxia-inducible factor-1alpha is associated with risk of aggressive behavior and tumor angiogenesis in gastrointestinal stromal tumor[END_REF].

Introduction

The integration of environmental factors influencing cell proliferation is the actual challenge of cell cycle modeling [START_REF] Csiksz-Nagy | Computational systems biology of the cell cycle[END_REF]. Indeed, since the 1990's, many models were developed to describe the evolution of protein levels during the cycle [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF][START_REF] Csikasz-Nagy | Analysis of a generic model of eukaryotic cell-cycle regulation[END_REF][START_REF] Barberis | Cell size at s phase 29 initiation: an emergent property of the g1/s network[END_REF][START_REF] Yang | Variations in cyclin d1 levels through the cell cycle determine the proliferative fate of a cell[END_REF]. However, few of them take into account external parameters, such as the temperature, the mechanical properties of the substrate, or hypoxia. This last factor is particularly interesting to study. In physiological conditions, hypoxia can occur if the blood does not bring enough oxygen to the cells. If hypoxia is too severe or too long, the cell enters into apoptosis [START_REF] Harris | Hypoxia-a key regulatory factor in tumour growth[END_REF]. In pathological conditions, cancer cells undergo a chronic hypoxia [START_REF] Kunz | Molecular responses to hypoxia in tumor cells[END_REF][START_REF] Lesart | On the importance of the submicrovascular network in a computational model of tumour growth[END_REF]. This hypoxia induces more aggressive, metastatic and resistant tumors [START_REF] Harris | Hypoxia-a key regulatory factor in tumour growth[END_REF][START_REF] Hernandez-Luna | Hypoxia inducible factor-1 induces chemoresistance phenotype in nonhodgkin lymphoma cell line via up-regulation of bcl-(xl)[END_REF]. In particular, if hypoxia can induce apoptosis in normal proliferating cells, tumor cells resist to apoptosis [START_REF] Harris | Hypoxia-a key regulatory factor in tumour growth[END_REF]. One mechanism that can explain this resistance is the entry into a quiescent state, where the cell stops its division cycle [START_REF] Royds | Response of tumour cells to hypoxia: role of p53 and nfkb[END_REF][START_REF] Alarcon | Quiescence: a mechanism for escaping the effects of drug on cell populations[END_REF]. This quiescent state is also a mean for the cancer cell to escape from the effects of chemotherapy [START_REF] Alarcon | Quiescence: a mechanism for escaping the effects of drug on cell populations[END_REF][START_REF] Letonqueze | Microrna-mediated posttranscriptional mechanisms of gene expression in proliferating and quiescent cancer cells[END_REF][START_REF] Powathil | Modelling the effects of cell-cycle heterogeneity on the response of a solid tumour to chemotherapy: biological insights from a hybrid multiscale cellular automaton model[END_REF]. Therefore, it participates to the aggressiveness of the tumors. Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF] proposed a simple model explaining why cancer cells can enter into quiescence under hypoxic conditions whereas normal cells follow their cycle. The authors introduce a difference between normal proliferating and cancer cells through the deregulation of the expression of a protein of the cell-cycle (p27) during tumorigenesis. However, this deregulation event is not a general feature of cancer. Indeed, this phenomenon was not observed in the majority of tumor cells [START_REF] Leshem | Phosphorylation of prb is required for hgfinduced muscle cell proliferation and is p27kip1-dependent[END_REF][START_REF] Park | Lovastatin-induced e2f-1 modulation and its effect on prostate cancer cell death[END_REF]. Besides, it is not clear that p27 is necessary to induce hypoxia-induced cell cycle arrest [START_REF] Box | Cell cycle kinase inhibitor expression and hypoxia-induced cell cycle arrest in human cancer cell lines[END_REF][START_REF] Green | p21(cip1) and p27(kip1) regulate cell cycle reentry after hypoxic stress but are not necessary for hypoxia-induced arrest[END_REF].

The HiF-1α factor is the central protein involved in the intracellular signaling pathway of hypoxia [START_REF] Semenza | Hypoxia-inducible factor 1 (hif-1) pathway[END_REF]. It is a transcription factor, which enhances the expression of numerous genes. These genes enable the cell to adapt to the environmental conditions (angiogenesis, arrest of aerobic metabolism), or to enter into apoptosis [START_REF] Goda | Hypoxia-inducible factor 1alpha is essential for cell cycle arrest during hypoxia[END_REF][START_REF] Papandreou | Hif-1 mediates adaptation to hypoxia by actively downregulating mitochondrial oxygen consumption[END_REF].

HiF-1 is an heterodimeric protein, constituted of two sub-units: HiF-1β, which is constitutively expressed, and HiF-1α, which is the sensor of hypoxia. Indeed, in order to activate hypoxia genes, HiF-1α has to be in a reduced form. When the level of oxygen is sufficient, an enzyme called HiF-1α prolylhydroxylase is active and converts HiF-1α into an hydroxylated form. This oxygenated form of HiF-1α is rapidly degraded by the proteasom pathway [START_REF] Harris | Hypoxia-a key regulatory factor in tumour growth[END_REF][START_REF] Salceda | Hypoxia-inducible factor 1alpha (hif-1alpha) protein is rapidly degraded by the ubiquitin-proteasome system under normoxic conditions. its stabilization by hypoxia depends on redox-induced changes[END_REF][START_REF] Powathil | Modelling the effects of cell-cycle heterogeneity on the response of a solid tumour to chemotherapy: biological insights from a hybrid multiscale cellular automaton model[END_REF].

The ways of actions of HiF-1α are numerous and complex [START_REF] Papandreou | Cellular reaction to hypoxia: sensing and responding to an adverse environment[END_REF]. A review of the huge molecular biology literature dealing with the effects of HiF-1α on the cell cycle enables us to retain several types of actions. First, it is clear that HiF-1α indirectly downregulates cyclin E activity, and this inhibition is the reason why HiF-1α causes slowing down or arrest of cell cycle [START_REF] Hammer | Hypoxic 32 suppression of the cell cycle gene cdc25a in tumor cells[END_REF][START_REF] Goda | Hif-1 in cell cycle regulation, apoptosis, and tumor progression[END_REF][START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF].

The origins of this action on cyclin E activity remain poorly defined [START_REF] Hammer | Hypoxic 32 suppression of the cell cycle gene cdc25a in tumor cells[END_REF], even if some potential pathways are known. Notably, the upregulation of cyclins inhibitors, such as p21 and p27, are reported [START_REF] Gardner | Hypoxia inhibits g1/s transition through regulation of p27 expression[END_REF][START_REF] Hammer | Hypoxic 32 suppression of the cell cycle gene cdc25a in tumor cells[END_REF][START_REF] Goda | Hypoxia-inducible factor 1alpha is essential for cell cycle arrest during hypoxia[END_REF]. However, some authors showed that the action of HiF-1α on p27 is not so clear since the expression of p27 under hypoxia may be independent of HiF-1α [START_REF] Box | Cell cycle kinase inhibitor expression and hypoxia-induced cell cycle arrest in human cancer cell lines[END_REF][START_REF] Carmeliet | Role of hif-1alpha in hypoxia-mediated apoptosis, cell proliferation and tumour angiogenesis[END_REF].

The second important effect of HiF-1α is the interrelation between this factor and cyclin D [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF][START_REF] Frei | Drosophila cyclin d/cdk4 requires hif-1 prolyl hydroxylase to drive cell growth[END_REF]. Wen et al. [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF] studied the effect of HiF stabilization on cyclin D level. They first found that after 24h at 0.2%, the mean cylin D concentration in the whole cell population decreases with 50% compared to the normoxic condition. To confirm the implication of HiF-1, its activity was impaired by DN-HiF overexpression, which induces an increase of the cyclin D level. There is also an action of cyclin D on HiF-1α, due to the activation of the HiF-1 prolyl-hydroxylase activity [START_REF] Frei | Drosophila cyclin d/cdk4 requires hif-1 prolyl hydroxylase to drive cell growth[END_REF]. Besides, the activation of cyclin D expression by HiF-2 underlines the rivalry in hypoxic tumor growth and progression between HiF-1 and HiF-2. Here, we choose to focus on the influence of cyclin D inhibition by HiF-1α, in order to show how this simple relationship can generate a serie of interesting results in agreement with data from the literature. Besides, the inhibition of the cyclins by HiF-1α under hypoxia is very well documented. As we previously said, hypoxia-dependent inhibition of cyclin E is considered as a cause of hypoxiadependent cell-cycle arrest. Second, the increase of unphosphorylated versus phosphorylated form of the Retinoblastoma protein [START_REF] Goda | Hif-1 in cell cycle regulation, apoptosis, and tumor progression[END_REF][START_REF] Goda | Hypoxia-inducible factor 1alpha is essential for cell cycle arrest during hypoxia[END_REF] under hypoxia is also a good evidence for the downregulation of cyclin D by HiF.

The increase of HiF-1α activity during tumorigenesis is well documented, since it seems to be a very common feature of cancers [START_REF] Zhong | Overexpression of hypoxia-inducible factor 1alpha in common human cancers and their metastases[END_REF][START_REF] Ke | Hypoxia-inducible factor-1 (hif-1)[END_REF][START_REF] Brown | Exploiting tumour hypoxia in cancer treatment[END_REF][START_REF] Rankin | The role of hypoxia-inducible factors in tumorigenesis[END_REF]. As a consequence, this factor became a new therapeutic target [START_REF] Semenza | Targeting hif-1 for cancer therapy[END_REF][START_REF] Semenza | Hypoxia-inducible factors: mediators of cancer progression and targets for cancer therapy[END_REF]. In some cases, its high level is simply due to the chronic hypoxia undergone by the tumor cells. In other cases, genetic mutations induce an over-stabilization of HiF-1. It is the case for the renal clear carcinoma [START_REF] Dai | Up-regulation of hypoxia inducible factor-1 by cobalt chloride correlates with proliferation and apoptosis in pc-2 cells[END_REF]. As a consequence, for the same oxygen pressure, HiF-1 will have a higher level in cancer cells than in normal cells [START_REF] Filippi | Csb protein is (a direct target of hif-1 and) a critical mediator of the hypoxic response[END_REF].

In a general way, a link was found between the aggressiveness of cancers and HiF-1α activity. This role of HiF in tumorigenesis is firstly due to the induction of anti-apoptotic and pro-angiogenic genes. Second, it mediates the entrance into quiescence of proliferating cells, which induces a hypoxiadependent chemoresistance.

The aim of this work is to describe an example of a simple mechanism of HiF-1α-dependent entrance into quiescence under hypoxia. We make a link between the level of HiF-1α and the ability of the cell to enter into quiescence under hypoxia. In agreement with the literature, our model assumes the regulation of HiF-1α stability by the oxygen pressure [START_REF] Winston | Growth factor regulation of cyclin d1 mrna expression through protein synthesis-dependent and -independent mechanisms[END_REF][START_REF] Salceda | Hypoxia-inducible factor 1alpha (hif-1alpha) protein is rapidly degraded by the ubiquitin-proteasome system under normoxic conditions. its stabilization by hypoxia depends on redox-induced changes[END_REF][START_REF] Maxwell | The pvhl-hif-1 system. a key mediator of oxygen homeostasis[END_REF]. We built a mathematical model of the G1/S transition in hypoxic conditions, which explicitly integrates the HiF-1α pathway. We focused on the relationships between HiF-1α and the cyclins, in order to propose a simple and biologically accurate mechanism of hypoxia-induced quiescence.

Model

Hypotheses

Our model is primary based on the models described in Alarcon et al.

[1], Tyson and Novak [START_REF] Tyson | Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions[END_REF], Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. In those papers, the G1/S transition is modeled by a biological switch between a cyclin and an inhibitor complex. They considered that this inhibitor was APC/cdh1. However, this role of APC/cdh1 complex is not well admitted; classicaly, the complex addressing cyclin E for degradation is known to be the SCF complex [START_REF] Van Drogen | Ubiquitylation of cyclin e requires the sequential function of scf complexes containing distinct hcdc4 isoforms[END_REF]. As a consequence, it is more accurate to consider that the inhibitor of cyclin E is SCF. The switch that we aim to study is represented in Figure 1. The G1 phase finishes when this molecular switch occurs. In our model, we differentiate the cyclin E and the cyclin D. The cyclin E is involved in the molecular switch where it inhibits and is inhibited by the SCF complex. Whereas Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] consider that two cyclins drive the G1/S transition (cyclin E and A), we consider here just one cyclin we named cyclin E. The concentration in active SCF complex is controlled by an evolution equation similar to that given by Tyson and Novak [START_REF] Tyson | Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions[END_REF] for the APC/cdh1 complex.

The cyclin D phosphorylates the Rb protein, which releases the transcription factor E2F. As Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF], we supposed that at each time, E2F Rb is in an equilibrium relationship with E2F and Rb because the complexation/decomplexation of phopshorylated Rb and E2F are supposed to be very fast compared to the evolution of the cyclin concentration. Besides, free E2F can be in an unphosphorylated (active) form, or in a phosphorylated (inactive) form. A dynamic equilibrium ensures the transition from phosphorylated to unphosphorylated forms. The total concentration of E2F (free phosphorylated/unphosphorylated, and complexed with Rb) is supposed to be constant, as in [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. Active E2F factor promotes the synthesis of the cyclin E, as assumed for the cyclin A by Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. Besides, to take into account the cell growth in the progression through the G1 phase, we assume that the mass increases the cyclin E concentration. This assumption is taken from Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF] and Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. The oxygen pressure is taken into account with the variable P , which represents the percentage of oxygen among all the other gases in the cell environment. The effect of the oxygen pressure is modeled by the concentration of HiF-1α. This concentration increases when the oxygen pressure decreases, as well described in the literature [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF][START_REF] Zhou | Assessment of hypoxia inducible factor levels in cancer cell lines upon hypoxic induction using a novel reporter construct[END_REF]. Jiang et al. [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF] show that in vitro cultured cells submitted to hypoxia exhibit an exponential relationship between the HiF-1 level and the oxygen pressure [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF]. The experimental results from Zhou et al. [START_REF] Zhou | Assessment of hypoxia inducible factor levels in cancer cell lines upon hypoxic induction using a novel reporter construct[END_REF] confirm that an exponential law is a good model for the hypoxiainduced HiF-1 stabilization. These authors also showed that this exponential law strongly depends on the cancer cell lines. We will discuss this point later.

HiF-1α decreases the cyclin D activity by inhibiting its synthesis, as observed in experiments [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF]. The evolution of the cyclin D level during the cycle, the influence of HiF-1α on its level, and the reaction of Rb phosphorylation are calibrated with data from the literature. We do not take into account the role of the p27 protein, contrary to Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF]. As we said, the effects of hypoxia on this protein remain unclear; besides, we want to focus on the influence of cyclin D/HiF interactions which promote the entrance into quiescence. Figure 2 gives a sketch of the molecular network.

Biological description

The variables used in the model are the cyclin D (cycD), the cyclin E (cycE), the Retinoblastoma protein (in its unphosphorylated form) (Rb), HiF-1α (H), the SCF complex (SCF ), the mass of the cell, and the different chemical forms of the E2F transcription factor. E2F exists in phosphorylated and unphosphorylated (E2F ) form, and free or linked to Rb (E2F Rb ). Thus, free E2F can be phosphorylated or unphosphorylated. The free and unphosphorylated form (E2F A ) is the active form. A chemical equilibrium drives the transition between the phosphorylated and unphosphorylated forms of E2F. We named E2F tot the total of all the chemical forms of E2F. As Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF], we assume this total to be constant. The cyclin D liberates E2F by phosphorylating Rb (that is to say by decreasing [Rb]). E2F A has a crucial role, since it promotes the synthesis of the cyclin E. Cell growth is the second factor influencing the cyclin E synthesis, trough the variable m.

HiF-1 inhibits the cyclin D synthesis. As a consequence, hypoxia induces a decrease of E2F A maximal level, which decreases the synthesis rate of the cyclin E. The core of the model is the molecular switch between the cyclin E and its inhibitor (SCF). Following the same idea as Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] and Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF], this switch is driven by the mutual inhibition of SCF and the cyclin E.

Mathematical formulation

The model can be written with a system of coupled ordinary differential equations (ODEs).

The cyclin D is synthesized with rate a 1 and naturally degraded with rate a 2 . We assume that it is linearly decreasing with [H] [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF]:

d[cycD] dt = a 1 -a 3 [H] -a 2 [cycD] (1)
The synthesis of cyclin E is induced by the cell growth and E2F A , according to Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] and Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF]. Its degradation is stimulated by the SCF complex:

d[cycE] dt = b 1 m[E2F A ] -b 2 [cycE] -b 3 [cycE][SCF ] (2)
The SCF complex obeys the same equation as the APC/cdh1 complex of the previous models [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF][START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF][START_REF] Tyson | Regulation of the eukaryotic cell cycle: molecular antagonism, hysteresis, and irreversible transitions[END_REF]:

d[SCF ] dt = e 1 (1 -[SCF ]) J 1 + 1 -[SCF ] -e 2 [SCF ][cycE] J 2 + [SCF ] (3)
As in Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF], Rb (unphosphorylated Rb) is synthetized with a d 2 rate, and naturally degraded with a d 3 rate. The cyclin D phosphorylates Rb. As a consequence, [Rb] decreases:

d[Rb] dt = d 2 -(d 3 + d 1 [cycD])[Rb] (4)
The growth of the cell mass follows the same law as in Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF],

with a growth rate α and a maximal size m 0 :

dm dt = αm(1 - m m 0 ) (5)
The unphosphorylated form of E2F follows a dynamic equilibrium described by Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]:

d[E2F ] dt = g 1 ([E2F ] tot -[E2F ]) (6)
The equilibrium relationship between the oxygen pressure and the HiF-1α concentration has been derived from experiments [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF]:

[H] = H 0 e β 1 (1-P ) (7)
Following Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF], we assume that at each time, the active E2F is the fraction of unphosphorylated

( [E2F ] [E2F ]tot ) free ([E2F ] tot -[E2F Rb ) E2F factor: [E2F A ] = ([E2F ] tot -[E2F Rb ])[E2F ] [E2F ] tot (8)
The equilibrium between E2F Rb , E2F and Rb is given by the equation from Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]:

[E2F Rb ] = 2[E2F ] tot [Rb] [E2F ] tot + [Rb] + ([E2F ] tot + [Rb]) 2 -4[E2F ] tot [Rb] (9)
As in previous models, the core of the system is the molecular switch between a cyclin and SCF. Mathematically, we consider that the cell passes through the G1 phase when [cycE] > [cycE] threshold and [SCF ] < [SCF ] threshold .

We chose [SCF ] thresold = 0.004 and [cycE] thresold = 0.15. When these values are reached, the cell divides.

Inhibition of the cyclin D by HiF-1α. Eq.( 1) and ( 7) describe the relationship between the cyclin D, HiF-1α and the oxygen pressure P in the cell environment. This model considers the effects of a long term (almost 48h) hypoxia.

Thus, the level of HiF is stabilized. The cyclin D is synthesized with rate a 1 , and degraded with rate a 2 . Its synthesis is inhibited by HiF-1α, which concentration depends on oxygen pressure following the expression given in eq. ( 7). The parameter a 3 represents the action of HiF on the synthesis of cyclin D. The parameter β 1 is determined from Jiang et al. [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF] (see Appendices A). The variable P stands for the oxygen pressure in the cell environment, expressed as a normalized percentage of all the gases. As suggested by Jiang et al. [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF], we consider that the reference pressure (P =1) is 6% O 2 . All the pressures will be normalized with this reference pressure. The parameters a 1 , a 2 and a 3 are estimated with data coming from Wen et al. [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF] and Yang et al. [START_REF] Yang | Variations in cyclin d1 levels through the cell cycle determine the proliferative fate of a cell[END_REF] (see Appendix A).

Control of E2F concentration. We denote [E2F ] tot the total concentration of E2F, E2F RB the complex E2F/Rb, E2F the unphosphorylated form of E2F and E2F A the free and unphosphorylated (active) form. We describe the equilibrium between these different variables from Eqs.( 6),( 8),( As in Novak and Tyson [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF], the dynamics of the inhibitor depends on the cyclin E, and on the transition between an active and an inactive form. These dependencies are expressed through Michaelis-Menten terms with enzymatic constants J 1 and J 2 respectively. The parameter e 1 promotes the transition from inactive to active forms, and e 2 the transition from active to inactive forms.

Parameters. The parameters we used are given in Table 1. Their values were estimated from previous models [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF][START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] and experimental data of molecular biology taken from the literature (see Appendix A). We will show that we can reduce the variability of all the parameters to the variability of e 2 and α, which control the dynamics of the system. These parameters can be easily adjusted to biological data.

Results

We simulated the G1/S transition by solving the system of equations Eqs.( 1)-( 9) for different values of P . We used a standard fourth order Runge-Kutta method, implemented in Matlab. The values of the parameters used in our first simulations are given in Table 1. We chose the undetermined parameters in order to have a G1 phase duration of about 400 minutes. This duration was chosen based on the work of [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. The rise of cyclins D and A levels, which was interpreted by Novak and Tyson, and by ourselves, as the signal for the G1/S transition, occurs about 6-7

hours after the beginning of the cycle. The initial conditions we used are given in Table 2. The initial values for the concentrations were chosen to have values of the same order of magnitude to that used by Alarcon et al.

[1] and Novak [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. We can notice that the initial concentrations do not influence the final state of the system (see Appendices B). We consider that the cell passes through the G1 phase when [cycE] > [cycE] threshold = 0.15 and [SCF ] < [SCF ] threshold = 0.04. These threshold conditions were chosen to capture as precisely as possible the moment of the switch between the cyclin E and the SCF complex. In the rest of the paper, the notation [X] ∞ stands for the value of [X] when the system reaches its steady state.

Hypoxia induces entrance into quiescence

Figure 3 shows simulations of the G1/S transition under three levels of hypoxia. When P diminishes, the molecular switch occurs with a delay which increases when hypoxia becomes more severe. For a certain value of P (here, 0.01), this switch does not occur any more. As in [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF], this inability for the cell to pass through the G1 phase can be interpreted as the entry into the quiescent state. P * denotes the minimal pressure which enables the cell to pass through the G1 phase.The parameters are given in Table 1, and P * = 0.03.

We also plotted the duration of the G1 phase as a function of P (Figure 4). As hypoxia becomes more intense, the G1 phase duration increases, in agreement with the literature [START_REF] Gardner | Hypoxia inhibits g1/s transition through regulation of p27 expression[END_REF]. For the lowest value of the pressure, we observe an asymptotic behavior, which finishes by the entrance into quiescence.

In order to mathematically characterize the final state of the cell (quiescent or proliferating), the stationary state solutions of the system, depending on P , are numerically calculated. The details of the calculation are presented in the Appendices B. Figure 5 presents the stable stationary states of the system as a function of P . When P < P * , the system exhibits three stable stationary states. With our initial conditions ([SCF ] high and [cycE] low), the equilibrium state of the system is reached with [SCF ] > [SCF ] threshold and

[cycE] < [cycE] threshold . This situation corresponds to the quiescent state.

When P > P * , the system has a unique stationary state. It is characterized by a high level of cycE and a low concentration of SCF , which leads the cell to achieve its cycle.

Influence of parameters variability

Previous observations enable us to suppose that the two main features of the cell response to hypoxia are the increased duration of the G1 phase and the capacity to enter into quiescence. The latter can be quantified by the value of P * (pressure for the bifurcation of the dynamics system). The core of the system is the molecular switch due to SCF and cycE interactions.

The variations of P have an influence on these dynamics by changing [E2F A ] and m. Since these two variables determine the synthesis rate of cycE, they affect the G1 phase duration and equilibrium values of SCF and cycE.

The equilibrium value of m never changes in our model. We call this maximal mass value m 0 , and we have m ∞ =m 0 . We can simulate its influence on the dynamics of the system by modifying b 1 , which buffers the variations of m 0 in the production rate of cycE. Therefore, we do not consider m 0 as a key parameter for the dynamics of the system. The variations of the growth rate, by modifications of the parameter α, just have an influence on the cell cycle duration.

Final concentrations of the SCF and cycE directly depend on [E2F

A ] ∞ . We can express [E2F A ] ∞ as a function of [cycD] ∞ , d 1 , d 2 (see Appendix B).
[cycD] ∞ depends on the parameters of Eq.( 1) and ( 7) (a 1 , a 2 , a 3 , β 1 ). As a consequence, the variability of all these parameters is totally buffered by the variability of b 1 , which controls the production rate of cycE.

As it ensures the coupling between SCF and cycE, the parameter e 2 also controls the dynamics of the molecular switch. If the coupling is weak, it is more difficult for the system to generate a molecular switch. As a consequence, we can also consider the parameter e 2 as a key factor influencing the behavior of the system.

Finally, the whole variability of the adjustment parameters can be simulated by the variability of the parameters e 2 , b 1 and α.

Figure 6 shows that the (b 1 , e 2 ) couple determines the value of P * . When e 2 decreases, the coupling between SCF and cycE becomes weak, and it is necessary to increase the hypoxic level to observe quiescence. As a consequence, P * decreases. In parallel, when b 1 decreases, the production rate of cycE at the steady state is less important. Therefore, the cell enters more easily into quiescence, and P * increases. We note that, for a given value of b 1 , the variations of e 2 are sufficient to simulate a large panel of values for P * . As a consequence, in order to adjust the model to a given value of this limit pressure, we can set b 1 and consider e 2 as an adjustment parameter. Finally, the parameters controlling the dynamics of the system are e 2 (coupling between cyclin E and SCF) and α (growth rate).

Discussion

This paper proposes a mathematical model of the G1/S transition in hypoxic conditions. Our model was calibrated with data from the literature, and the results from other mathematical models of the G1 phase transition [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF][START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. The core of this model is a molecular switch between a cyclin and an inhibitor (SCF). The passage through the G1 phase is determined by the possibility to realize this switch. If the switch does not occur, we consider that the cell is trapped into a quiescent state. Hypoxia was simulated by introducing HiF-1α. We modeled the inhibition between cyclin D and HiF-1α, in agreement with data from the literature [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF]. The oxygen pressure determines, in this way, the final concentration of cyclins D during the G1 phase. This level of cyclins D determines the maximum value of the active E2F concentration, which controls the synthesis rate of cyclins E. If this rate is not high enough, the critical concentration for the cyclin E is not reached and the molecular switch cannot occur. The cell enters into a quiescent state.

As a consequence, if the oxygen pressure is low, the cyclin D will not reach a sufficient level to ensure the minimal amount of cyclins E the cell needs to achieve the G1 phase. We could define a limit pressure P * , which is the minimum oxygen pressure whereby the cell can pass through the G1 phase.

Mathematically, it corresponds to the bifurcation point of our system.

Thus, this model illustrates a simple mechanism of hypoxia-mediated slowing down of the cell cycle, and entrance into quiescence. It focuses on the influence of the interactions between HiF-1α and the cyclins on the G1 phase achievement. Notably, we show how the inhibition of cyclin D by HiF-1α can affect the cell cycle. By inhibiting the cyclin D activity, HiF promotes the entrance into quiescence. Thus, it provides a resistance to apoptosis for the proliferative cells. This result gives an explanation to the results of Wen et al. [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF], which have shown that the HiF-mediated chemoresistance appearing in tumor cells is linked to the inhibition of the cyclin D. Indeed, they observed that the inhibition of cyclin D expression by siRNA induces an increase of the resistance to chemically-induced apoptosis. They concluded than the negative regulation of the cyclin D induced by HiF can, at least partially, explain the chemoresistance which is associated to cancerous HiF-1 overstabilization. Thus, our model shows how a certain form of chemoresistance can be associated to HiF-1 in cancer cells. In the case of a similar response to hypoxia between normal and cancer cells (non-mutated hypoxia pathway), the model highlights a possible mechanism for the chronic hypoxia-induced entry into quiescence.

In the case of a genetic over-stabilization of HiF-1, which is the case for various type of cancers [START_REF] Dai | Up-regulation of hypoxia inducible factor-1 by cobalt chloride correlates with proliferation and apoptosis in pc-2 cells[END_REF][START_REF] Brdos | Negative and positive regulation of hif-1: a complex network[END_REF], the increased ability of cancer cells to enter into a quiescent state is well described by our model. Indeed, this case can be modeled by multiplying [H 0 ] by an ǫ > 1 parameter to artificially increase the HiF-1 level. As a consequence, for a same oxygen pressure, HiF-1 concentration will be more important in cancer cells than in normal cell.

Thus, for a given hypoxic condition, a cancer cell can enter into a quiescent (resistant) state whereas the normal cells continue their cycle (more sensitive to apoptosis). Our model thus captures the differential behavior of cancer and normal cell. We can notice than this ǫ parameter could depend on the cell line, since Zhou et al. [START_REF] Zhou | Assessment of hypoxia inducible factor levels in cancer cell lines upon hypoxic induction using a novel reporter construct[END_REF] showed that different cancer cell lines does not have the same normoxic level of HiF-1.

Compared to the model from Alarcon et al. [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF], our work describes in a very general way the influence of HiF-1 on hypoxia-induced quiescence of normal and cancer cells. Whereas the previous model was based on an arguable action of hypoxia on p27, we constructed a biologically more accurate system of ODEs. This new model gives a possible explanation to the link between HiF-1, the cyclins and hypoxia-induced chemoresistance. Besides, our model can easily be completed to study the influence of a genetic overstabilization of HiF-1 by considering an ǫ parameter, which can vary with the considered cell line.

The α parameter, which drives the growth speed of the cell, was chosen to obtain a value of 400min for the G1 phase duration. It is important to notice that the results of our model are not influenced by the choice of the normoxic G1 phase duration. Since it does not influence the final state of the mathematical system, this α parameter can be modified to adjust the model to various G1 phase durations. It is biologically correct, because the duration of the G1 phase in the absence of an environmental stress is known to be, at least for a part, driven by the growth rate of the cell [START_REF] Barberis | Cell size at s phase 29 initiation: an emergent property of the g1/s network[END_REF].

Besides, the influence of HiF-1 on the cell metabolism was not taken into account in this model. It is known that the glycolitic switch described by [START_REF] Gatenby | Why do cancers have high aerobic glycolysis?[END_REF] [START_REF] Gatenby | Why do cancers have high aerobic glycolysis?[END_REF] is linked to the activity of HiF-1. Notably, it was shown that HiF-1 stimulates the synthesis of glycolitic enzymes, such as Pyruvate Kinase [START_REF] Papandreou | Hif-1 mediates adaptation to hypoxia by actively downregulating mitochondrial oxygen consumption[END_REF]. In parallel, Icreverzi et al. [START_REF] Icreverzi | Drosophila cyclin d/cdk4 regulates mitochondrial biogenesis and aging and sensitizes animals to hypoxic stress[END_REF] observed that the cyclin D enhances the cell growth through the stimulation of mitobiogenesis. By inhibiting the synthesis of cyclins D, HiF-1 makes the mitobiogenesis to decrease, which promotes the transition to a glycolytic phenotype. This switch allows cancer cells to compensate low aerobic metabolism by an increased glucose consumption. This can be a way for the cell to escape the quiescence state. It could be interesting to add to our model a term describing the energetic production, and to link the cell growth to this available energy.

We also have to notice that our model only studies the effect of a given oxygen pressure on a single cell. In the reality of tumor growth, there is a gradient of oxygen which varies in space and time. To study the effects of such a dynamical hypoxic condition, it would be necessary to consider some series of cycles, by resetting the values of the variables to their initial values after the G1/S transition is reached. For each new cycle, the oxygen level can change, and the cell response may also change.

Besides, it is possible to extend this model to consider a population of interacting cells. These interactions can be, for instance, through the diffusion of growth factors which stimulate cell proliferation by enhancing the synthesis of cyclins D [START_REF] Winston | Growth factor regulation of cyclin d1 mrna expression through protein synthesis-dependent and -independent mechanisms[END_REF]. We can also consider the differential mechanical interactions between the tumor cells, and between tumor cells and the stroma. These mechanical interactions are known to influence cell division [START_REF] Provenzano | Mechanical signaling through the cytoskeleton regulates cell proliferation by coordinated focal adhesion and rho gtpase signaling[END_REF].

We can notice here that the kinetics of HiF-1 accumulation is not considered in this model. We study the idealized case of a constant hypoxia, and we consider the effect of the equilibrium value of HiF-1. However, these kinetics were measured ( [START_REF] Moroz | Real-time of hif-1alpha stabilization and degradation[END_REF][START_REF] Kamat | Mutant p53 facilitates pro-angiogenic, hyperproliferative phenotype in response to chronic relative hypoxia[END_REF]). HiF-1 shows a peak after 6h, and then decreases to stabilize at an equilibrium level after 24-48h. This equilibrium level depends on the intensity of hypoxia, and on the cellular type (normal and cancer cells). If we consider the more complex case of a dynamical hy-poxia, it should be necessary to take into account these kinetics of HiF-1 accumulation.

We note that the exact relationships between HiF-1 and the cyclin D is not totally elucidated. Many reports show that HiF-1 inhibits the cyclin D [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF][START_REF] Frei | Drosophila cyclin d/cdk4 requires hif-1 prolyl hydroxylase to drive cell growth[END_REF] and the cyclin E [START_REF] Hammer | Hypoxic 32 suppression of the cell cycle gene cdc25a in tumor cells[END_REF][START_REF] Goda | Hif-1 in cell cycle regulation, apoptosis, and tumor progression[END_REF][START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF] leading to the cell cycle arrest. However, an other report [START_REF] Baba | Loss of von hippel-lindau protein causes cell density dependent deregulation of cyclind1 expression through hypoxia-inducible factor[END_REF] shows that HiF-1 can also stimulate the induction of cyclins D, leading to cell proliferation. In our work, we use the first biological hypothesis, and we show that it can explain the cell cycle arrest phenomenon that is commonly observed.

a 3 H 0 0.0085 [50] β 1 2.5 [27] b 1 0.018 b 2 0.5 [1] b 3 1 [1] J 1 0.04 [1] J 2 0.04 [1] d 1 0.2 [36] d 2 0.1 [1] e 2 14
e 1 1 α 0.005 g 1 0.016 [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] [E2F] tot 1

[cycE] threshold 0.15 [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF] [SCF ] threshold 0.004 [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF]

1 .

 1 ously proposed by Novak and Tyson [36]. The first describes the transition, mediated by g 1 , from phosphorylated to unphosphorylated E2F. The others are equilibrium equations modeling the relationship between the different chemical species involving E2F. Eq.(4) describes the evolution of the unphosphorylated Retinoblastoma protein concentration. [Rb] is synthesized with coefficient d 2 , naturally phosphorylated in the cell environment with coefficient d 3 , and phosphorylated by the cyclin D with coefficient d Molecular switch. Eqs.(2),(3) describe the biological switch which forms the core of the model for the G1/S transition. Besides, we suppose that the synthesis of cyclins E, represented by the variable cycE in our model, depends on the mass m and on the transcription factor E2F in its active form (E2F A )[START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF]. The parameter b 1 regulates the speed of the chemical synthesis. In the other hand, cycE is naturally degraded with coefficient b 2 , and suffers SCFmediated degradation with coefficient b 3 . The aim is to simulate a switch similar to that obtained in these previous models of the G1/S transition.

Finally, we aim

  to calibrate the undetermined parameters with experimental data, by measuring the G1 phase duration under hypoxia and the value of P * . Since we showed that only two parameters (α and e 2 ) drive the dynamics of the system (G1 phase duration and entrance into a quiescence state), experimental data will enable us to set the values of these parameters for different cell lines. Therefore, this would provide a tool to compare different cancer cell lines on the basis of their ability to enter into a quiescence state. It could be interesting to build up a classification of different cancer cells using the values of the parameters of our model. By comparing it with a known classification made on the basis of their resistance and aggressiveness, we could obtain some information about the correlation between the ability to enter into quiescence under hypoxia and the aggressiveness.

  Figure 1:

Figure 3

 3 Figure 3:

Figure 4

 4 Figure 4:

Figure 1 :

 1 Figure 1: The biological switch between cyclin E and its inhibitor is necessary to pass through the G1/S restriction point. Schematic diagram showing the evolution of the concentration of SCF and cyc-CDK during the cell cycle. The vertical dotted line represents the restriction point: when the cell achieves the switch, it finishes the cycle. From [1].

Figure 2 :

 2 Figure 2: Model of the G1/S transition in hypoxic conditions: schematic representation of the molecular network considered in our model. We modeled the action of cyclin D, which phosphorylates Rb and liberates E2F. This transcription factor activates the synthesis of cyclin E, which is degraded by SCF complex. HiF-1α negatively regulates cyclin D synthesis. RbNP=Nonphosphorylated Rb.

Figure 3 :

 3 Figure 3: Evolution of [SCF ] (solid black line) and [cycE] (solid grey line). The simulations were made for adimensionalized values of P , with P =1(top), 0.08 (center), 0.01 (bottom). When the molecular switch occurs, the cell passes through the G1 phase. For P =0.01, the system does not generate a switch any more. The cell enters into a quiescent state.

Figure 4 :

 4 Figure 4: Influence of hypoxia on the G1 phase duration. Value of the G1 phase duration (in min) given by the model as a function of the parameter P (adimensionalized).

Figure 5 :

 5 Figure 5: Influence of P (adimensionalized) on the stable fixed points of the system. [SCF ] ∞ (top) and [cycE] ∞ (bottom) are plotted as functions of P * . The dark line represent the stationary states reach by the system withthe parameters we used. The point where we pass from three to one stable stationary state corresponds to P * (limit of the quiescent state).

Figure 6 :

 6 Figure 6: Influence of b 1 (synthesis rate of cycE) and e 2 (coupling between cycE and SCF ) on the value of the bifurcation point P * . The value of P * (color scale) is plotted as a function of b 1 and e 2 . The white zones correspond to non physiological values (the cell never or always finishes its cycle).

Table 1 :

 1 Parameters used in the model. The parameters we could not estimate are initially

	determined in order to reach a G1 phase duration of about 400 min for proliferative cells
	in normoxic conditions

Table 2 :

 2 Initial values taken for the simulations

Table legends Table 1 :

 legends1 Parameters used in the model. The parameters we could not estimate are initially determined in order to reach a G1 phase duration of about 400 min for proliferative cells in normoxic conditions.

Table 2 :

 2 Initial values taken for the simulations.

Acknowledgements

We are very grateful to Françoise Giroud for her advices which helped us to build and justify the hypotheses of the model.

Appendix A. Determination of the parameters

The parameter we did not adjust were taken from [START_REF] Alarcon | A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells[END_REF] and [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] models, or calculated from the molecular biology literature. This first part presents the biological data used to build the model.

Appendix A.1. HiF-1α and oxygen pressure

In [START_REF] Jiang | Hypoxia-inducible factor 1 levels vary exponentially over a physiologically relevant range of o2 tension[END_REF] the relationship between [H] (HiF-1 concentration) and P was experimentally determined. The authors have shown that the level of HiF-1α varies with oxygen pressure in an exponential form, when the oxygen pressure is below 6%. The experimental curve they present was fitted in order to build the dimensionless expression describing the relationship between [H] and P .

The normoxia (P =1) is reached for an oxygen pressure P =6%.

Appendix A.2. Hif-1α and cyclin D dynamic

In [START_REF] Stacey | Cyclin d1 serves as a cell cycle regulatory switch in actively proliferating cells[END_REF], the level of cyclin D is quantified during the cycle: the protein is multiplied by a factor taken between 2 and 5 during the G1 phase. We adjust the maximum level of cyclin D to reach a concentration five times higher than in the beginning of the G1 phase. We also used the dynamic of the cyclin D evolution during the G1 phase modeled by [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF], with a rapid increase of cyclin D concentration (maximum level reached in less than one hour). In [START_REF] Wen | Suppression of cyclin d1 by hypoxia-inducible factor-1 via direct mechanism inhibits the proliferation and 5-fluorouracil-induced apoptosis of a549 cells[END_REF] experiments were realized to quantify the inverse correlation between HiF-1α and cyclin D. Their paper shows that after 20 hours of hypoxia (oxygen pressure P =0.21%) in a cell constitutively expressing the cyclin D, its concentration reaches an equilibrium, and is divided by about 2 with respect to normoxia. Considering that our model describes cells that have endured several hours of hypoxia, and that the constitutive expression is equivalent to the induction at the beginning of the G1 phase, the parameter a 3 was adjusted to obtain a maximal value of [cycD] divided by 2 when P = 0.21%.

Appendix A.3. Phosphorylation of Rb

We used the model developed by [START_REF] Novak | A model for restriction point control of the mammalian cell cycle[END_REF] to adjust the parameters controlling the quantity of phosphorylated Rb: they obtained a phosphorylation of 50% of total Rb after 4 hours of cycle.

Appendix B. Fixed points

X ∞ denotes the steady state concentration of the variable X. By equaling to zero Eqs.( 1)-( 6), it is possible to calculate the equilibrium values of the variables used in the model. The calculation of cycD ∞ , Rb ∞ , E2F ∞ , m ∞ is straightforward:

By equaling to zero Eqs.( 2) and (3), we obtain an equation for cycE ∞ :

The equation B.2 is solved numerically, for different values of P and of the parameters α, b 1 , e 2 . The resolution gives us the number and the values of the fixed points; we also numerically determined the stability of these fixed points.


A Polyakov
email: andrey.polyakov@inria.fr

Discontinuous Lyapunov Functions for Nonasymptotic Stability Analysis

The paper presents mathematical tools required for finite-time stability analysis of discontinuous control systems using discontinuous Lyapunov functions. Elements of Filippov theory of differential equations with discontinuous right-hand sides and stability notions are briefly observed. Concepts of generalized derivatives and non-smooth Lyapunov functions are considered. The generalized Lyapunov theorems for stability analysis and convergence time estimation are presented and supported by examples.

INTRODUCTION

During whole history of control theory, a special interest of researchers was focused on systems with relay (switching) control elements [START_REF] Flugge-Lotz | Discontinuous automatic control[END_REF], [START_REF] Tsypkin | Relay Control Systems[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Boiko | Discontinuous control systems : frequencydomain analysis and design[END_REF]). Such control systems have found applications in many engineering areas. They are simple, effective, cheap and sometimes they have better dynamics than linear systems [START_REF] Tsypkin | Relay Control Systems[END_REF]). In practice both input and output of a system may be of a relay type. For example, automobile engine control systems sometimes use λ -sensor with almost relay output characteristics, i.e the only sign of a controllable output can be measured [START_REF] Choi | Robust throttle control of automotive engines[END_REF]). In the same time, terristors can be considered as relay "actuators" for some power electronic systems [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF]).

Mathematical backgrounds for a rigorous study of relay systems were presented in the beginning of 1960s by the celebrated Filippov's theory of differential equations with discontinuous right-hand sides, see [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]. Following this theory, discontinuous differential equations have to be extended to differential inclusions by means of a special procedure. It is worth to stress that Filippov's theory was severely criticized by many authors, see [START_REF] Neimark | Note on A. Filippov paper[END_REF], [START_REF] Aizerman | Foundations of theory of discontinuous systems i[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], since it does not describe adequately some discontinuous and relay models. That is why, extensions and specifications of this theory appear rather frequently, see, for example, [START_REF] Bartolini | Variable structure systems nonlinear in the control law[END_REF], [START_REF] Heemels | Input-to-state stability and interconnections of discontinuous dynamical systems[END_REF].

Analysis of relay systems is frequently related to a specific property, which is called finite-time stability [START_REF] Roxin | On finite stability in control systems[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF], [START_REF] Orlov | Finite time stability and robust control synthesis of uncertain switched systems[END_REF], [START_REF] Moulay | Finite-time stability conditions for non-autonomous continuous systems[END_REF]). Indeed, the simplest example of a finite-time stable system is the relay system of the form: ẋ =sign[x], x ∈ R, x(0) = x 0 . Any solution of this system reaches the origin in a finite time T (x 0 ) = |x 0 | and remains there for all later time instants. Sometimes, this conceptually very simple property is hard to proof theoretically. From a practical point of view, it is also important to estimate a time of stabilization (settling time). Both these problems can be tackled by Lyapunov Function Method [START_REF] Lyapunov | The general problem of the stability of motion[END_REF], [START_REF] Barbashin | Lyapunov Functions[END_REF], [START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF]). However, designing a finite-time Lyapunov function of a rather simple form is a difficult problem for many systems. For concrete relay systems appropriate Lyapunov functions can be nonsmooth (Polyakov and Poznyak [2009b,a]), non-Lipschitz [START_REF] Moreno | Strict Lyapunov functions for the super-twisting algorithm[END_REF]) and even discontinuous [START_REF] Polyakov | Unified Lyapunov function for a finite-time stability analysis of relay second-order sliding mode control systems[END_REF]. Some problems of asymptotic and finitetime stability analysis using Lipschitz continuous Lyapunov functions are studied in, for example, [START_REF] Shevitz | Lyapunov stability theory of nonsmooth systems[END_REF], [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF], [START_REF] Bacciotti | Nonpathological Lyapunov functions and discontinuous caratheodory systems[END_REF], [START_REF] Moulay | Finite-time stability conditions for non-autonomous continuous systems[END_REF]. The lower semicontinuous Lyapunov functions for asymptotic stability analysis are discussed in [START_REF] Chellaboina | Generalized Lyapunov and invariant set theorems for nonlinear dynamical systems[END_REF], [START_REF] Clarke | Qualitative properties of trajectories of control systems: A survey[END_REF], [START_REF] Zhu | Lower semicontinuous lyapunov functions and subdifferential calculus[END_REF] This paper elaborates the non-asymptotic (finite-time) stability analysis of differential inclusions based on discontinuous Lyapunov functions. The main tool to be applied for this purpose is the theory of contingent derivatives (see [START_REF] Zaremba | Sur une extension de la notion d'equation differentiel[END_REF], [START_REF] Natanson | Theory of functions of a real variable[END_REF], [START_REF] Mordukhovich | Variational Analysis and Generalized Differentiation I: Basic Theory[END_REF]).

The paper is organized as follows. The next section presents notations are used in the paper. Then the elements of the Filippov theory differential equations with discontinuous right-hand sides and differential inclusions are discussed. Next the evolution of the stability notion is observed. After then the supporting results about contingent derivatives are studied in order to present a discontinuous Lyapunov function method in the last section.

The following notations are used in the paper R is the set of real numbers and

R = R ∪ {-∞} ∪ {+∞}, R + = {x ∈ R : x > 0} and R + = R + ∪ {+∞}; I denotes one of the following intervals: [a, b], (a, b), [a, b) or (a, b],
where a, b ∈ R, a < b; the set consisting of elements x 1 , x 2 , ..., x n is denoted by {x 1 , x 2 , ..., x n }; the set of all subsets of a set M ⊆ R n is denoted by 2 M ; the inner product of x, y ∈ R n is denoted by x, y and x =

x, x ; the sign function is defined by

sign σ [ρ] = 1 if ρ > 0, -1 if ρ < 0, σ if ρ = 0, (1)
where σ ∈ R : -1 ≤ σ ≤ 1. If σ = 0 we use the notation sign [ρ]; the set-valued modification of the sign function is given by sign

[ρ] = {1} if ρ > 0, {-1} if ρ < 0, [-1, 1] if ρ = 0;
(2) the geometric sum of two sets is denoted by " +", i.e.

M 1 +M 2 = x1∈M1,x2∈M2 {x 1 + x 2 }, (3)
where ; W n I is the set of vector-valued, componentwise locally absolutely continuous functions, which map I to R n .

M 1 ⊆ R n , M 2 ⊆ R n ; the

SYSTEM DESCRIPTION

Let us consider the system ẋ ∈ F (t, x), t ∈ R, (4) where x ∈ R n is a system state, F : R × R n → 2 R n is a setvalued mapping, which satisfies the standard assumptions: the set F (t, x) is nonempty, convex and compact for any

(t, x) ∈ R × R n ; the set-valued mapping F is upper semicontinuous in R × R n .
The conventional approach to the regularization of the differential equation with discontinuous right-hand side ẋ

= f (t, x), t ∈ R, x ∈ R n , f ∈ C -1
(5) is its extension to the differential inclusion. For instance, the Filippov procedure gives the following construction for the right-hand side of the differential inclusion (4)

F (t, x) = δ>0 µ(N )=0 co f (t, {x} +B(δ)\N ), (6)
where the intersections are taken over all sets N ⊂ R n of measure zero (µ(N ) = 0) and all δ > 0. The Filippov regularization provides the differential inclusions, which satisfy the standard assumptions. Definition 1. [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF], page 50). An absolutely continuous function x : I → R n defined on some interval or segment I is called a solution of (5) if it satisfies the differential inclusion (4), ( 6) almost everywhere on I. Theorem 1. [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF], page 77). Let a set-valued function F : G → 2 R n be defined and upper semi-continuous at each point of the set

G = {(t, x) ∈ R n+1 : |t - t 0 | ≤ a and x -x 0 ≤ b}, where a, b ∈ R + , t 0 ∈ R, x 0 ∈ R n . Let F (t,
x) be nonempty, compact and convex for (t, x) ∈ G.

If there exists K > 0 such that y < K for all y ∈ F (t, x) and all (t, x) ∈ G then there exists at least one absolutely continuous function x : R → R n defined at least on the segment [t 0α, t 0 + α], α = min{a, b/K}, such that x(t 0 ) = x 0 and the inclusion ẋ(t) ∈ F (t, x(t)) holds almost everywhere on [t 0α, t 0 + α].

Alternative regularization procedures for discontinuous control systems can be found, for example, in [START_REF] Aizerman | Foundations of theory of discontinuous systems i[END_REF], [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Heemels | Input-to-state stability and interconnections of discontinuous dynamical systems[END_REF].

STABILITY AND CONVERGENCE RATE

Consider the differential inclusion (4) for t > t 0 with an initial condition x(t 0 ) = x 0 , (7) where x 0 ∈ R n is given.

Cauchy problem ( 4), ( 7) may have non-unique solution for a given t 0 ∈ R and a given x 0 ∈ R n . Let us denote the set of all solutions of Cauchy problem (4), ( 7) by Φ(t 0 , x 0 ) and a solution of (4), ( 7) by x(t, t 0 , x 0 ) ∈ Φ(t 0 , x 0 ). Nonuniqueness of solutions implies two types of stability for differential inclusions (4): weak stability(a property holds for a solution) and strong stability (a property holds for all solutions) (see, for example, [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF]). Weak stability usually is not enough for robust control purposes. This section observes only strong stability properties of the system (4). All conditions presented in definitions below are assumed to be held for all solutions x(t, t 0 , x 0 ) ∈ Φ(t 0 , x 0 ).

Lyapunov and asymptotic stability

The concept of stability introduced in the famous thesis of A.M. Lyapunov (see, its recent translation in [START_REF] Lyapunov | The general problem of the stability of motion[END_REF]) is the central notion of the modern stability theory. It considers some nominal motion x * (t, t 0 , x 0 ) of a dynamic system and studies small perturbations of the initial condition x 0 . If they imply small deviations of perturbed motions from x * (t, t 0 , x 0 ) then the nominal motion is called stable. We study different stability forms of the zero solution (or, equivalently, the origin) of the system (4), since making the change of variables y = xx * it is nothing to transform any problem of stability analysis for some nontrivial solution x * (t, t * , x * 0 ) to the same problem for the zero solution. Assume that 0 ∈ F (t, 0), ∀t ∈ R, where F is defined by (4). Definition 2. (Lyapunov stability). The origin of the system ( 4) is said to be Lyapunov stable if for ∀ε ∈ R + and ∀t 0 ∈ R there exists δ = δ(ε, t 0 ) ∈ R + such that : 1) any solution x(t, t 0 , x 0 ) (with x 0 ∈ B(δ)) of Cauchy problem (4), ( 7) exists for all t > t 0 ; 2)

x(t, t 0 , x 0 ) ∈ B(ε) for t > t 0 if x 0 ∈ B(δ).
If the function δ does not depend on t 0 then the origin is called uniformly Lyapunov stable. For instance, if F (t, x) is independent of t (time-invariant case) and the zero solution of (4) is Lyapunov stable, then it is uniformly Lyapunov stable. The origin, which does not satisfy any condition from Definition 2, is called unstable. Proposition 1. If the origin of (4) is Lyapunov stable then x(t) = 0 is the unique solution of (4), ( 7) with x 0 = 0. Definition 3. (Asymptotic attractivity). The origin of the system (4) is said to be asymptotically attractive if for ∀t 0 ∈ R there exists a set U (t 0 ) ⊆ R n : 0 ∈ int(U (t 0 )) : 1) any solution x(t, t 0 , x 0 ) (with x 0 ∈ U(t 0 )) of Cauchy problem (4), ( 7) exists for t > t 0 ; 2) lim t→+∞ x(t, t 0 , x 0 ) = 0. The set U (t 0 ) is called attraction domain. Definition 4. (Asymptotic stability). The origin of the system ( 4) is said to be asymptotically stable if it is Lyapunov stable and asymptotically attractive.

If U (t 0 ) = R n then the asymptotically stable (attractive) origin of the system (4) is called globally asymptotically stable (attractive).

Requirement of Lyapunov stability is very important in Definition 4, since even global asymptotic attractivity does not imply Lyapunov stability (see, for example, [START_REF] Vinograd | The inadequacy of the method of characteristic exponents for the study of nonlinear differential equations[END_REF], page 433 or [START_REF] Hahn | Stability of Motion[END_REF], page 191).

Finite-time Stability

Introduce the functional

T 0 : W n [t0,+∞) → R + ∪ {0} by the following formula T 0 (y(•)) = inf τ ≥t0:y(τ )=0 τ. If y(τ ) = 0 for all t ∈ [t 0 , +∞) then T 0 (y(•)) = +∞.
The settling-time function of the system (4) let us define by

T (t 0 , x 0 ) = sup x(t,t0,x0)∈Φ(t0,x0) T 0 (x(t, t 0 , x 0 )) -t 0 , (8)
where Φ(t 0 , x 0 ) is the set of all solutions of the Cauchy problem ( 4), ( 7). Definition 5. (Finite-time attractivity). The origin of the system ( 4) is said to be finite-time attractive if for ∀t 0 ∈ R there exists a set V(t 0 ) ⊆ R n : 0 ∈ int(V(t 0 )) : 1) any solution x(t, t 0 , x 0 ) (with x 0 ∈ V(t 0 )) of Cauchy problem ( 4), ( 7) exists for all t > t 0 ; 2) T (t 0 , x 0 ) < +∞ for x 0 ∈ V(t 0 ) and for t 0 ∈ R. The set V(t 0 ) is called finite-time attraction domain.

It is worth to stress that the presented finite-time attractivity property, introduced originally in Bhat and Bernstein [2000], is "weaker" than asymptotic attractivity. However, it is important for many control applications. For example, antimissile control problem has to be studied only on a finite interval of time, since there is nothing to control after missile explosion. In practice, Lyapunov stability is additionally required in order to guarantee a robustness of a control system. Definition 6. (Finite-time stability [START_REF] Roxin | On finite stability in control systems[END_REF])). The origin of the system (4) is said to be finite-time stable if it is Lyapunov stable and finite-time attractive.

If V(t 0 ) = R n then the origin of ( 4) is called globally finitetime stable. Example 1.

Consider the system ẋ ∈ - 2 √ π • sign[x] +{|2tx|}, t > t 0 , x ∈ R, ( 9
)
where t 0 ∈ R. It can be shown that the origin of this system is finite-time attractive with an attraction domain

V(t 0 ) = B e t 2 0 (1 -erf(|t 0 |))
, where erf(z) = 2 √ π z 0 e -τ 2 dτ , z ∈ R is the so-called Gauss error function. Moreover, the origin of the considered system is Lyapunov stable (for ∀ε > 0 and for ∀t 0 ∈ R we can select δ = δ(t 0 ) = min ε, e t 2 0 (1erf(|t 0 |)) ), so it is finite-time stable. In particular, for t 0 > 0 the settling-time function has the form T (t 0 , x 0 ) = erf -1 |x 0 |e -t 2 0 + erf(t 0 )t 0 , where erf -1 (•) denotes the inverse function to erf(•).

The proposition 1 implies the following claim. Proposition 2. [START_REF] Bhat | Finite time stability of continuous autonomous systems[END_REF]). If the origin of the system (4) is finite-time stable then it is asymptotically stable and x(t, t 0 , x 0 ) = 0 for t > t 0 + T 0 (t 0 , x 0 ).

GENERALIZED DERIVATIVES

The celebrated Second Lyapunov Method is founded on the so-called energetic approach to stability analysis. It considers any positive definite function as an possible energetic characteristic ("energy") of a dynamic system and studies evolution of this "energy" in time. If a dynamic system has an energetic function, which is decreasing (strictly decreasing) along any trajectory of the system, then this system has a stability property and the corresponding energetic function is called Lyapunov function.

For example, to analyze asymptotic stability of the origin of the system ẋ

= f (t, x), f ∈ C(R n+1 ), t ∈ R + , x ∈ R n (
10) it is sufficient to find a continuous positive definite function V (•) such that for any solution x(t) of the system (10) the function V (x(t)) is decreasing and tending to zero for t → +∞. The existence of such function guarantees asymptotic stability of the origin of the system (10) due to Zubov's theorem (see, [START_REF] Poznyak | Advanced Mathematical Tools for Automatic Control Engineers[END_REF]).

If the function V (x) is continuously differentiable then the required monotonicity property can be rewritten in the form of the classical condition [START_REF] Lyapunov | The general problem of the stability of motion[END_REF]

: V (x(t)) = ∇ T V (x)f (t, x) < 0.
(11) The inequality ( 11) is very usable, since it does not require knowing the solutions of (10) in order to check the asymptotic stability. From the practical point of view, it is important to represent monotonicity conditions in the form of differential or algebraic inequalities like (11).

Finite-time stability analysis of relay systems is frequently based on non-smooth or even discontinuous Lyapunov functions [START_REF] Roxin | On finite stability in control systems[END_REF], Polyakov and Poznyak [2009b], [START_REF] Moreno | Strict Lyapunov functions for the super-twisting algorithm[END_REF], which require consideration of generalized derivatives and generalized gradients in order to verify stability conditions. This section presents all necessary backgrounds for the corresponding non-smooth analysis.

Derivative Numbers and Monotonicity

Let I be one of the following intervals: Recall that the function ϕ : R → R is called decreasing on

I iff ∀t 1 , t 2 ∈ I : t 1 ≤ t 2 ⇒ ϕ(t 1 ) ≥ ϕ(t 2 ).
Let K be a set of all sequences of real numbers converging to zero, i.e. {h n } ∈ K ⇔ h n → 0, h n = 0. Let a real-valued function ϕ : R → R be defined on I. Definition 7. [START_REF] Natanson | Theory of functions of a real variable[END_REF], page 207). A number D {hn} ϕ(t)= lim n→+∞ ϕ(t+hn)-ϕ(t) hn

, {h n } ∈ K : t + h n ∈ I is called derivative number of the function ϕ at a point t ∈ I, if finite or infinite limit exists.

The set of all derivative numbers of the function ϕ at the point t ∈ I is called contingent derivative:

D K ϕ(t) = {hn}∈K {D {hn} ϕ(t)} ⊆ R.
A contingent derivative of a vector-valued function ϕ : R → R n can be defined in the same way. If a function ϕ(t) is differentiable at a point t ∈ I then D K ϕ(t) = { φ(t)}. Lemma 1. [START_REF] Natanson | Theory of functions of a real variable[END_REF], page 208). If a function ϕ : R → R is defined on I then 1) a set D K ϕ(t) ⊆ R is nonempty for any t ∈ I;

2) for any t ∈ I and for any sequence {h n } ∈ K : t + {h n } ∈ I there exists a subsequence {h n ′ } ⊆ {h n } such that finite or infinite derivative number D {h n ′ } ϕ(t) exists.

Lemma 1 remains true for a vector-valued function ϕ : R → R n . Inequalities y < 0, y ≤ 0, y > 0, y ≥ 0 for y ∈ R n are understood in a componentwise sense. If for ∀y ∈ D K ϕ(t) we have y < 0 then we write D K ϕ(t) < 0.

Other ordering relations ≤, >, ≥ for contingent derivatives are interpreted analogously.

The contingent derivative also helps to prove monotonicity of a non-differentiable function.

Lemma 2. [START_REF] Natanson | Theory of functions of a real variable[END_REF], page 266). If a function ϕ : R → R is defined on I and the inequality D K ϕ(t) ≤ 0 holds for all t ∈ I, then ϕ(t) is decreasing function on I and differentiable almost everywhere on I.

Lemma 2 require neither the continuity of the function ϕ(t) nor the finiteness of its derivative numbers. It gives a background for the discontinuous Lyapunov function method.

Example 2. The function ϕ(t) = -tsign σ [t] has a negative contingent derivative for all t ∈ R and for any σ ∈ [-1, 1], where the function sign σ is defined by (1). Indeed,

D K ϕ(t) = {-1} for t = 0, D K ϕ(0) = {-∞} if σ ∈ (-1, 1) and D K ϕ(0) = {-∞, -1} if σ ∈ {-1, 1}.
The next lemma simplifies the monotonicity analysis of nonnegative functions. Lemma 3. If 1) the function ϕ : R → R is nonnegative on I; 2) the inequality D K ϕ(t) ≤ 0 holds for t ∈ I : ϕ(t) = 0; 3) the function ϕ(t) is continuous at any t ∈ I : ϕ(t) = 0; then ϕ(t) is decreasing function on I and differentiable almost everywhere on I.

The lemma of Wazewski given below presents the relation between solutions of a differential inclusion (4) and its contingent derivatives. Lemma 4. [START_REF] Filippov | Differential Equations with Discontinuous Righthand Sides[END_REF], page 70). Let a set-valued function F : R n+1 → 2 R n be defined, upper-semicontinuous on a closed nonempty set Ω ∈ R n+1 and the set F (t, x) be nonempty, compact and convex for all (t, x) ∈ Ω. Let an absolutely continuous function x : R → R n be defined on I and (t,

x(t)) ∈ Ω if t ∈ I. Then ẋ(t) ∈ F (t, x(t)) almost everywhere on I ⇔ D K x(t) ⊆ F (t, x(t))
everywhere on I.

Generalized directional derivatives

Stability analysis based on Lyapunov functions requires calculation of derivatives of positive definite functions along trajectories of a dynamic system. If Lyapunov function is non-differentiable, a concept of generalized directional derivatives (see, for example, [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]) can be used for this analysis. This survey introduces generalized directional derivatives by analogy with contingent derivatives for scalar functions.

Let M(d) be a set of all sequences of real vectors converging to

d ∈ R n , i.e. {v n } ∈ M(d) ⇔ v n → d, v n ∈ R n .
Let a function V : R n → R be defined on an open nonempty set Ω ⊆ R n and d ∈ R n . Definition 8. A number

D {hn},{vn} V (x, d) = lim n→+∞ V (x + h n v n ) -V (x) h n , {h n } ∈ K, {v n } ∈ M(d) : x + h n v n ∈ Ω is called directional derivative number of the function V (x) at the point x ∈ Ω on the direction d ∈ R n , if finite or infinite limit exists.
The set of all directional derivative numbers of the function V (x) at the point x ∈ Ω on the direction

d ∈ R n is called directional contingent derivative: D K,M(d) V (x) = {hn}∈K,{vn}∈M(d) {D {hn},{vn} V (x, d)}.
Similarly to Lemma 1, if x ∈ Ω then the set D K,M(d) V (x) is nonempty for any function V defined on an open nonempty set Ω ⊆ R n and any d ∈ R n . The next lemma gives a chain rule for the introduced contingent derivative. Lemma 5. Let a function V : R n → R be defined on an open nonempty set Ω ⊆ R n and a function x : R → R n be defined on I, such that x(t) ∈ Ω if t ∈ I and the contingent derivative

D K x(t) ⊆ R n is bounded for all t ∈ I. Then the inclusion D K V (x(t)) ⊆ d∈D K x(t) D K,M(d) V (x) holds ∀t ∈ I.
This lemma together with Lemmas 4 and 3 imply the corollary, which is useful for Lyapunov analysis. Corollary 1. Let a set-valued function F : R n+1 → 2 R n be defined and upper-semicontinuous on I × Ω and a set F (t, x) be nonempty, compact and convex for any (t, x) ∈ I × Ω, where Ω ⊆ R n is an open nonempty set. Let x(t, t 0 , x 0 ) be an arbitrary solution of Cauchy problem (4), (7) defined on [t 0 , t 0 + α), where t 0 ∈ I, x 0 ∈ Ω and α ∈ R + . Let a function V : R n → R be nonnegative on Ω. If the inequality D F (t,x) V (x) ≤ 0 holds for every t ∈ I and every x ∈ Ω : V (x) = 0 then the function of time V (x(t, t 0 , x 0 )) is decreasing on [t 0 , t 0 + α), where

D F (t,x) V (x) = d∈F (t,x) D K,M(d) V (x).
(12)

DISCONTINUOUS LYAPUNOV FUNCTION METHOD FOR FINITE-TIME STABILITY ANALYSIS

Lyapunov function method is a very effective tool for analysis and design of both linear and nonlinear control systems [START_REF] Bacciotti | Lyapunov Functions and Stability in Control Theory[END_REF]. Initially, the method was presented for "unrated" (Lyapunov and asymptotic) stability analysis [START_REF] Lyapunov | The general problem of the stability of motion[END_REF]. A development of control theory had required to study a convergence rate together with a stability properties of a control system. This section presents the discontinuous Lyapunov function method for finite-time stability analysis. Let us refund initially the Zubov's theorems about Lyapunov and asymptotic stability using the contingent derivative approach.

Analysis of Lyapunov and asymptotic stability

Recall that a continuous function W : R n → R defined on R n is said to be positive definite iff W (0) = 0 and W (x) > 0 for x ∈ R n \{0}. Definition 9. A function V : R n → R is said to be proper on an open nonempty set Ω ⊆ R n : 0 ∈ int(Ω) iff 1) it is defined on Ω and continuous at the origin; 2) there exists a continuous positive definite function V : R n → R such that V (x) ≤ V (x) for x ∈ Ω.

A positive definite function W : R n → R is called radially unbounded if W (x) → +∞ for x → +∞. Definition 10. A function V : R n → R is said to be globally proper iff it is proper on R n and the positive definite function V : R n → R is radially unbounded.

If V is continuous on Ω, then V (x) = V (x) for x ∈ Ω and Definition 10 corresponds to the usual notion of proper positive definite function (see [START_REF] Clarke | Asymptotic stability and smooth Lyapunov functions[END_REF]).

For a given number r ∈ R and a given positive definite function W : R n → R defined on Ω let us introduce the set Π(W, r) = {x ∈ Ω : W (x) < r}, which is called the level set of the function W .

Theorems on Lyapunov and asymptotic stability given below are obtained by a combination of Zubov's theorems (see, for example, [START_REF] Poznyak | Advanced Mathematical Tools for Automatic Control Engineers[END_REF], pages 566-568) with Corollary 1. Theorem 2. Let a function V : R n → R be proper on an open nonempty set Ω ⊆ R n : 0 ∈ int(Ω) and D F (t,x) V (x) ≤ 0 for t ∈ R and x ∈ Ω\{0}.

(13) Then the origin of the system (4) is Lyapunov stable.

Asymptotic stability requires analysis of an attraction set. Lyapunov function approach may provide an estimate of this set. Theorem 3. Let a function V : R n → R be proper on an open nonempty set Ω ⊆ R n : 0 ∈ int(Ω), a function W : R n → R be a continuous positive definite and D F (t,x) V (x) ≤ -W (x) for t ∈ R and x ∈ Ω\{0}. Then the origin of the system (4) is asymptotically stable with an attraction domain If V is globally proper and Ω = R n then the origin of the system (4) is globally asymptotically stable (U = R n ).

U = Π(V, λ(h)) ∩ B(h), (14)
If V is locally bounded on Ω then the origin of the system (4) is uniformly asymptotically stable.

The presented theorems shows that discontinuous and non-Lipschitz Lyapunov functions can also be used for stability analysis using the theory of contingent derivatives.

The following important theorem declares that a smooth Lyapunov function always exists for a time-invariant asymptotically stable differential inclusion (4). Theorem 4. [START_REF] Clarke | Asymptotic stability and smooth Lyapunov functions[END_REF], Theorem 1.2). Let a set-valued function F : R n → R n be defined and uppersemicontinuous in R n . Let F (x) be nonempty, compact and convex for any x ∈ R n . If the origin of the system ẋ ∈ F (x) is globally uniformly asymptotically stable iff there exists a globally proper function

V (•) ∈ C ∞ (R n ) and a function W (•) ∈ C ∞ (R n ) : W (x) > 0 for x = 0 such that max y∈F (x) ∇ T V (x)y ≤ -W (x), x ∈ R n \{0}.
However, the practice shows that designing of a Lyapunov function for nonlinear and/or discontinuous system is a nontrivial problem even for a two dimensional case.

Frequently, in order to analyze stability of a sliding mode control system it is simpler to design initially a non-smooth Lyapunov function (see, for example, [START_REF] Utkin | Sliding Mode Control in Electro-Mechanical Systems[END_REF], [START_REF] Polyakov | Unified Lyapunov function for a finite-time stability analysis of relay second-order sliding mode control systems[END_REF], [START_REF] Moreno | Strict Lyapunov functions for the super-twisting algorithm[END_REF]) and after then to use some technique for its smoothing, see, for example, [START_REF] Clarke | Asymptotic stability and smooth Lyapunov functions[END_REF].

Lyapunov analysis of finite-time stability

Analysis of finite-time stability using the Lyapunov function method allows us to estimate of a settling time a priori. The proof of the next theorem follows the ideas introduced in [START_REF] Roxin | On finite stability in control systems[END_REF], [START_REF] Polyakov | Unified Lyapunov function for a finite-time stability analysis of relay second-order sliding mode control systems[END_REF].

Theorem 5. Let a function V : R n → R be proper on an open nonempty set Ω ⊆ R n : 0 ∈ int(Ω) and D F (t,x) V (x) ≤ -q for t ∈ R and x ∈ Ω\{0}, (15) for some q ∈ R + . Then the origin of the system (4) is finite-time stable with an attraction domain U defined by ( 14) and T (x 0 ) ≤ V (x 0 )/q for x 0 ∈ U, where T (•) is a settling-time function.

If a function V is globally proper on Ω = R n then the inequality (15) implies global finite-time stability of the system (4).

Obviously, a settling-time function of finite-time stable system is a Lyapunov function in a generalized sense. Theorem 5 operates with a very large class Lyapunov functions. However, its conditions are still rather conservative. For example, the settling-time function from the next example does not satisfy Theorem 5. However, its simple discontinuous modification transforms the corresponding settling time function to a Lyapunov function in the sense of Theorem 5. Example 3. Consider the so-called twisting system [START_REF] Levant | Sliding order and sliding accuracy in sliding mode control[END_REF])

ẋ1 ẋ2 ∈ F (x 1 , x 2 ) = y -(2sign[x 1 ] +sign[x 2 ]) , (16)
which is finite-time stable with the settling-time function [START_REF] Polyakov | Unified Lyapunov function for a finite-time stability analysis of relay second-order sliding mode control systems[END_REF]):

T tw (x 1 , x 2 ) = p |x 1 | + x 2 2 2(2 + sign[x 1 x 2 ]) + |x 2 | sign[x 1 x 2 ] 2 + sign[x 1 x 2 ] ,
where p = 4 √ 2/(3 -√ 3). The function T tw is globally proper, Lipschtz continuous outside the origin and continuously differentiable for x 1 x 2 = 0 D F T tw = ∂Ttw ∂x1 x 2 -∂Ttw ∂x2 (2 sign[x 1 ]+ sign[x 2 ])= -1. (17) In the same time, T tw (x 1 , x 2 ) does not satisfy (15), since D F T tw (0, x 2 ) R + = ∅ for x 1 = 0. Applying generalized Clark's gradient (see, for example, [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]) does not solves this problem. The paper presents the discontinuous Lyapunov function method for non-asymptotic stability analysis of differential inclusions. The example of discontinuous Lyapunov function is presented for twisting control system in order to demonstrate some advances of the method.

  [a, b], (a, b), [a, b) or (a, b], where a, b ∈ R, a < b.

  where λ(h) = inf x∈R n : x =h V (x) and h = sup r∈R+:B(r)⊆Ω r.

However

  Fig. 1. A trajectory and the discontinuous Lyapunov function of twisting system 6. CONCLUSIONS

then the function T tw (x 1 , x 2 ) is discontinuous on the line x 1 = 0 for x 2 = 0, but it is globally proper and satisfies the conditions of Theorem 5. Indeed, the equality (17) still holds and D F T tw (0, x 2 ) = {-∞} for x 2 = 0, D F T tw (x 1 , 0) ≤ -1/3 for x 1 = 0. The Fig.3depicts a trajectory and the discontinuous Lyapunov function for twisting system.