domain: sdv math

text:
Swarup Satya

Aurélien Samal
email: samal@cs.uni-bonn.de

Dima Naldi

Andreas Grigoriev

Nathalie Weber
email: weber@cs.uni-bonn.de

Ovidiu Théret

Radulescu

Samal

Aurélien Naldi
email: aurelien.naldi@gmail.com

Dima Grigoriev
email: dmitry.grigoryev@math.univ-lille1.fr

Andreas Weber

Nathalie Théret
email: nathalie.theret@univ-rennes1.fr

Ovidiu Radulescu
email: ovidiu.radulescu@univ-montp2.fr

Geometric analysis of pathways dynamics: application to versatility of TGF-β receptors

Keywords: Tropical Geometry, Cancer Systems Biology, Finite State Automaton, Metastability

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Networks of biochemical reactions are used in computational biology as models of signalling, metabolism, and gene regulation. For various applications it is important to understand how the dynamics of these models depend on internal parameters, initial data and environment variables. Traditionally, the dynamics of biochemical networks is studied in the framework of chemical kinetics that can be either deterministic (ordinary differential equations) or stochastic (continuous time Markov processes). In order to cope with qualitative data, boolean [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF][START_REF] Thomas | Boolean formalization of genetic control circuits[END_REF] or multi-valued networks [START_REF] Thomas | Regulatory networks seen as asynchronous automata: a logical description[END_REF][START_REF] Naldi | Decision diagrams for the representation and analysis of logical models of genetic networks[END_REF] are used instead of continuous models. For this reason, many efforts were focused on coarse graining dynamical networks described by ordinary differential equations (ODE) to boolean networks [START_REF] Davidich | The transition from differential equations to boolean networks: a case study in simplifying a regulatory network model[END_REF]. However, in spite of their advantages, dynamical properties of large Boolean or multi-valued networks are still difficult to study. The difficulty originates from the number of possible states, which for multi-valued networks with m levels (Boolean networks correspond to m = 2) is m n . Although there are efficient methods for computing attractors of Boolean networks (algorithms based on binary decision diagrams or on satisfiability solvers can handle synchronous networks with hundreds of variables [START_REF] Dubrova | A sat-based algorithm for finding attractors in synchronous boolean networks[END_REF]), more intricate dynamical properties like the metastable regimes discussed in this paper, ask for comprehensive analysis of the state transition graph which is hard to perform and moreover for analysis of the hierarchy of time scales which is even harder for Boolean and multi-valued networks. The coarse graining method proposed in this paper leads to a drastic reduction of the number of states. This offers unprecedented possibilities for qualitative analysis of the dynamics.

In this paper we propose a new method for model analysis that uses coarse grained descriptions of continuous dynamics as discrete automata defined on finite states. These states will not be obtained by discretization of network variables, but by identification of a finite number of collective modes describing possible coordinated activity of several variables.

For large networks with ordinary differential equations dynamics and multiple timescales it is reasonable to consider the following property: a typical trajectory consists in a succession of qualitatively different slow segments separated by faster transitions. The slow segments, generally called metastable states or regimes, can be of several types such as attractive invariant manifolds [START_REF] Gorban | Invariant manifolds for physical and chemical kinetics[END_REF], Milnor attractors [START_REF] Rabinovich | Dynamical principles in neuroscience[END_REF] or saddle connections [START_REF] Rabinovich | Information flow dynamics in the brain[END_REF]. The notion of metastability generalizes the notion of attractor. Like in the case of attractors, distant parts of the system can have coordinated activity for metastability. The dynamical states of large networks can be represented as points in a high dimensional space, called phase space. In this representation each coordinate represents the concentration of a molecular species. Coordinated activity means that many of the species concentrations are correlated, which can be geometrically represented by the proximity to a lower dimension hypersurface in the phase space. A system remains in the proximity of an attractor after entering its basin of attraction, but can leave a metastable regime after a relatively long time (much longer than the time needed for transitions between two different regimes). Fig. 1 summarizes this geometrical picture. The term crazy-quilt was coined to describe such a patchy landscape of multiscale networks dynamics [START_REF] Gorban | Dynamic and static limitation in reaction networks, revisited[END_REF].

Fig. 1 here This phenomenon, called itinerancy, received particular attention in neuroscience [START_REF] Tsuda | Chaotic itinerancy as a dynamical basis of hermeneutics in brain and mind[END_REF]. Itinerant behaviour is shown by mathematical models with stable heteroclinic sequences (defined as open chains of saddle fixed points connected by one-dimensional separatrices) and was also observed in transient activity of antennal lobe neurons involved in insect olfaction or in the activity of high vocal centers controlling songbird patterns [START_REF] Rabinovich | Dynamical principles in neuroscience[END_REF]. We believe that similar phenomena occur in molecular biology for chemical reaction networks. A well studied example sustaining this picture is the biochemical and gene expression dynamics guiding the orderly progression of the cell cycle. The main feature of this dynamics is the sequential activation of cyclin dependent kinase/cyclin complexes. More than 30 years since cyclins were discovered the main cell cycle control events are now well understood and it is agreed that each of them involve the collective action of several regulator molecules. In addition, studies of periodic gene expression in synchronized cell division of yeast show the existence of waves of coordinate expression corresponding to different cell cycle phases or transitions [START_REF] Rustici | Periodic gene expression program of the fission yeast cell cycle[END_REF]. Furthermore, mathematical models of the cell cycle machinery (currently more than 150 published models [START_REF] Weis | A datadriven, mathematical model of mammalian cell cycle regulation[END_REF]) illustrate the stage dependent coordination of biochemical variables. As an example, the structure of steady state branches of the Wee1-Cdc25 module in yeast lead [START_REF] Tyson | The dynamics of cell cycle regulation[END_REF] to consider that the exit from mitosis is a collective phenomenon that can be described as a saddle-node bifurcation. Our analysis of such models also showed that branches and bifurcations of states occur naturally in the context of cell cycle modelling [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF]. In a more general context, geometric analysis of single-cell expression data from human and mouse tissues showed that gene expression is structured in clusters but also in continua of states within polyhedra whose vertices can be understood as specialized key tasks [START_REF] Korem | Geometry of the gene expression space of individual cells[END_REF]. These findings were interpreted in terms of multi-objective optimization solutions [START_REF] Korem | Geometry of the gene expression space of individual cells[END_REF], but could also suggest transient behaviour between specialized states. The idea of associating cell lineage commitment to collective behaviour of gene networks variables was used in many other contexts including cancer genomics where it was proposed that cancer cells are trapped in some abnormal attractors [START_REF] Huang | Cancer attractors: A systems view of tumors from a gene network dynamics and developmental perspective[END_REF].

In this paper we propose a method to compute metastable dynamical regimes and the transitions between such regimes for chemical reaction networks. This will provide a precise meaning to the "crazy-quilt" metaphor illustrated in Fig. 1. To this aim we will use tropical geometry methods. Tropical methods [START_REF] Litvinov | Maslov dequantization, idempotent and tropical mathematics: a brief introduction[END_REF][START_REF] Maclagan | Introduction to tropical geometry[END_REF], also known as idempotent or max-plus algebras owe their name to the fact that one of the pioneer of the field, Imre Simon, was Brazilian. These methods found numerous applications to computer science [START_REF] Simon | Recognizable sets with multiplicities in the tropical semiring[END_REF], physics [START_REF] Litvinov | Maslov dequantization, idempotent and tropical mathematics: a brief introduction[END_REF], railway traffic [START_REF] Chang | On deterministic traffic regulation and service guarantees: a systematic approach by filtering[END_REF], and statistics [START_REF] Pachter | Tropical geometry of statistical models[END_REF]. Recently we have applied these methods to model order reduction [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF][START_REF] Noel | Tropicalization and tropical equilibration of chemical reactions[END_REF]Radulescu et al., 2015b;Samal et al., 2015b). In these works we have used tropical methods to rank monomial terms into rate vectors according to their orders of magnitude and to identify lowest order, dominant terms. When there is only one dominant term or when the domi-nant terms have all the same sign, the dynamics is fast and the system tends rapidly towards a region in phase space where at least two dominant terms of opposite signs are equilibrated. We have called the latter situation tropical equilibration [START_REF] Noel | Tropicalization and tropical equilibration of chemical reactions[END_REF]Radulescu et al., 2015b;Samal et al., 2015b). In this paper, we use tropical equilibrations to identify metastable dynamic regimes of chemical reaction networks. We show that tropical equilibrations can be grouped into branches and describe the qualitative network dynamics as a sequence of transitions from one branch to another. The complexity of the qualitative dynamics depends on the number of branches. Therefore, we would like to know how this number depends on the number of chemical species. Although there are theoretical results suggesting that the number of branches should be small, these results are valid in the average in the probabilistic space of all the models. In order to test this property numerically we will compute the branches for a large collection of models of the Biomodels database [START_REF] Le Novere | BioModels database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF]. The idea to gather tropical equilibrations into branches was introduced in our CMSB 2015 conference paper (Radulescu et al., 2015a). As new work with respect to this paper, we propose the benchmarking of the method on the Biomodels database. We also investigate the biological significance of metastable states in the context of TGFβ signalling in cancer cell lines.

The structure of the paper is as follows. In the Theory and Methods section we introduce the branches of tropical equilibrations and discuss briefly how they can be calculated. We also propose an algorithm to learn Markov state models defined on branches of tropical equilibrations. The Results section has two parts. First, we apply the computation of branches to models in the Biomodels database. Then, we apply the method to a model of TGFβ signalling and show how the analysis can be used to interpret biological data.

Theory and Methods

Tropical equilibrations of chemical reactions networks with polynomial rate functions

In this section we introduce the main concepts relating geometry and dynamics.

We consider chemical reaction networks described by mass action kinetics. To describe such a reaction network we need the list of components, A = {A 1 , ... A n } and the list of reactions (the reaction mechanism):

i α j i A i k β j k A k , (1)
where j ∈ [1, r] is the reaction number, α j i and β j k are positive integers defining the reaction stoichiometry.

Then, the mass action kinetics is described by a system of differential equations:

dx i dt = r j=1 k j S ij x α j , 1 ≤ i ≤ n, (2)
where k j > 0 are kinetic constants, x i are variable concentrations, S ij = β j i -α j i are the entries of the stoichiometric matrix, α j = (α j 1 , . . . , α j n ) are multi-indices, and

x α j = x α j 1 1 . . . x α j n n .
For our reasonings, we can replace the exact values of parameters by their orders of magnitude that are supposed to be known. Usually, orders of magnitude are approximations of the parameters by integer powers of ten and serve for rough comparisons. Our definition of orders of magnitude is based on the equation k j = kj ε γ j , where ε is a small positive number. The exponents γ j are considered to be integer or rational. For instance, the approximation γ j = round(log(k j )/ log(ε)),

(3) produces integer exponents, whereas γ j = round(d log(k j )/ log(ε))/d produces rational exponents, where round stands for the closest integer (with half-integers rounded to even numbers) and d is a strictly positive integer. When ε = 1/10, our definition provides the usual decimal orders.

In this study, orders of magnitude of the kinetic parameters are supposed to be known. In contrast, species orders vary in time and have to be computed. To this aim, the species concentrations are first represented by orders of magnitude defined as a j = log(x j )/ log(ε).

(4) More precisely, one has x j = xj ε a j , where xj has zero order (unity). Because log(ε) < 0, Eq.( 4) means that species orders and concentrations are anticorrelated (large orders mean small concentrations and vice versa).

Then, network dynamics is described by a rescaled ODE system

dx i dt = j ε µ j (a)-a i kj S ij xα j , (5)
where µ j (a) = γ j + a, α j (6), and , stands for the dot product. The r.h.s. of each equation in ( 5) is a sum of multivariate monomials in the concentrations. The orders µ j indicate how large these monomials are in absolute value. A monomial of order µ j dominates another monomial of order µ j if µ j < µ j .

To set these ideas down let us use a simple chemical network example, the Michaelis-Menten kinetics:

S + E k 1 k -1 ES k 2 → P + E,
where S, ES, E, P represent the substrate, the enzyme-substrate complex, the enzyme and the product, respectively.

After using the two conservation laws E + ES = e 0 and S + ES + P = s 0 , we find ẋ1

= -k 1 x 1 (e 0 -x 2 ) + k -1 x 2 , ẋ2 = k 1 x 1 (e 0 -x 2 ) -(k -1 + k 2 )x 2 . ( 7
)
where x 1 , x 2 are the concentrations of S and ES respectively. Orders of variables and parameters are as follows

x i = xi ε a i , 1 ≤ i ≤ 2, k 1 = k1 ε γ 1 , k -1 = k-1 ε γ -1 , e 0 = ē0 ε γe .
The tropical equilibration problem consists in the equality of the orders of at least two monomials one positive and another negative in the differential equations of each species. This condition allows us to compute the concentration orders defined by (4). More precisely, we want to find a vector a such that min

j,S ij >0 (γ j + a, α j ) = min j,S ij <0 (γ j + a, α j ) (8)
The equation( 8) is related to the notion of tropical hypersurface. A tropical hypersurface is the set of vectors a ∈ R n such that the minimun min j,S ij =0 (γ j + a, α j ) is attained for at least two different indices j (with no sign conditions). Tropical prevarieties are finite intersections of tropical hypersurfaces. Therefore, our tropical equilibrations are subsets of tropical prevarieties [START_REF] Maclagan | Introduction to tropical geometry[END_REF]. The sign condition in (8) was imposed because species concentrations are real positive numbers. A sum of positive monomials can not be zero for positive real values of the variables.

The system (8) can be seen as a system of equations in min-plus algebra (also known as tropical semiring), where multiplication ⊗ is the real numbers addition x ⊗ y = x + y and the addition ⊕ is the minimum operation x ⊕ y = min(x, y). In order to find the solutions of such system we can explore combinatorially trees of solutions resulting from various choices of minimal terms and write down inequations for each situation. Because a set of inequations define a polyhedron, the set of tropical equilibration solutions forms a polyhedron in R n . As a matter of fact, computing tropical equilibrations from the orders of magnitude of the model parameters is a NP-complete problem [START_REF] Theobald | On the frontiers of polynomial computations in tropical geometry[END_REF] and brute force calculation by exploration of combinatorics has exponential complexity. However, methods based on the Newton polytope [START_REF] Samal | A tropical method based on newton polygon approach for algebraic analysis of biochemical reaction networks[END_REF] or constraint logic programming [START_REF] Soliman | A constraint solving approach to model reduction by tropical equilibration[END_REF] exploit the sparseness and redundancy of the system and reduce the combinatorics and the time to compute tropical solutions.

The tropical equilibration equations for the Michaelis-Menten example are obtained by equating minimal orders of positive monomials with minimal orders of negative monomials in (7):

γ 1 + γ e + a 1 = min(γ 1 + a 1 , γ -1 ) + a 2 , (9) γ 1 + γ e + a 1 = min(γ 1 + a 1 , min(γ -1 , γ 2 )) + a 2 . ( 10
)
Species timescales. The timescale of a variable x i is given by 1

x i dx i dt = 1 xi dx i dt whose order is ν i = min{µ j |S ij = 0} -a i . (11)
The order ν i indicates how fast is the variable

x i (if ν i < ν i then x i is faster than x i ) .
Partial tropical equilibrations. It is useful to extend the tropical equilibration problem to partial equilibrations, that means solving (8) only for a subset of species. This is justified by the fact that slow species do not need to be equilibrated. In order to have a self-consistent calculation we compute the species timescales by (11). A partial equilibration is consistent if ν i < ν for all non-equilibrated species i. ν > 0 is an arbitrarily chosen threshold indicating the timescale of interest.

Tropical equilibrations, slow invariant manifolds and metastable dynamic regimes. In dissipative systems, fast variables relax rapidly to some low dimensional attractive manifold called invariant manifold [START_REF] Gorban | Invariant manifolds for physical and chemical kinetics[END_REF] that carries the slow mode dynamics. A projection of dynamic equations onto this manifold provides the reduced dynamics [START_REF] Maas | Simplifying chemical kinetics: intrinsic lowdimensional manifolds in composition space[END_REF]. This simple picture can be complexified to cope with hierarchies of invariant manifolds and with phenomena such as transverse instability, excitability and itineracy. Firstly, the relaxation towards an attractor can have several stages, each with its own invariant manifold. During relaxation towards the attractor, invariant manifolds are usually embedded one into another (there is a decrease of dimensionality) [START_REF] Chiavazzo | Adaptive simplification of complex multiscale systems[END_REF]. Secondly, invariant manifolds can lose local stability, which allow the trajectories to perform large phase space excursions before returning in a different place on the same invariant manifold or on a different one [START_REF] Haller | Localized instability and attraction along invariant manifolds[END_REF]. We showed elsewhere that tropical equilibrations can be used to approximate invariant manifolds for systems of polynomial differential equations [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF][START_REF] Noel | Tropicalization and tropical equilibration of chemical reactions[END_REF]Radulescu et al., 2015b). Indeed, tropical equilibration are defined by the cancelling out of dominant forces acting on the system. The remaining weak non-compensated forces ensure the slow dynamics on the invariant manifold. Tropical equilibrations are thus different from steady states, in that there is a slow dynamics. In this paper we will use them as proxies for metastable dynamic regimes.

More precisely, let us assume that species timescales (defined by Eq.( 11)) satisfy the relation ν 1 ≤ ν 2 ≤ . . . ≤ ν n and that not all the timescales are the same, i.e. there is m < n such that ν m+1 -ν m > 0. Then, two groups of variables have separated timescales. The variables X r = (x 1 , x 2 , . . . , x m ) are fast (change significantly on timescales of order of magnitude ε -νm or shorter. The remaining variables X s = (x m+1 , x m+2 , . . . , x n ) are much slower (have little variation on timescales of order of magnitude ε -νm ). The metastable regime means that fast variables have reached quasi-steady state values on a low dimensional hypersurface of the phase space.

Branches of tropical equilibrations and connectivity graph. For each equation i, let us define

M i (a) = argmin j (µ j (a), S ij > 0) = argmin j (µ j (a), S ij < 0), (12)
in other words M i denotes the set of monomials having the same minimal order µ i . We call tropically truncated system the system obtained by pruning the system (5), i.e. by keeping only the dominating monomials.

dx i dt = ε µ i -a i ( j∈M i (a) kj ν ji xα j ), (13)
The tropical truncated system is uniquely determined by the index sets M i (a), therefore by the tropical equilibration a. Reciprocally, two tropical equilibrations can have the same index sets M i (a) and truncated systems. We say that two tropical equilibrations a 1 , a 2 are equivalent iff

M i (a 1 ) = M i (a 2 ), for all i. Equivalence classes of tropical equilibrations are called branches. A branch B with an index set M i is minimal if M i ⊂ M i for all i where M i is the index set of the branch B implies B = B or B = ∅.
Closures of equilibration branches are defined by a finite set of linear inequalities, which means that they are polyhedral complexes. Minimal branches correspond to maximal dimension faces of the polyhedral complex. The incidence relations between the maximal dimension faces (n -1 dimensional faces, where n is the number of variables) of the polyhedral complex define the connectivity graph. More precisely, minimal branches are the vertices of this graph. Two minimal branches are connected if the corresponding faces of the polyhedral complex share a n -2 dimensional face. In terms of index sets, two minimal branches with index sets M and M are connected if there is an index set M of an existing non-minimal branch such that M i ⊂ M i and M i ⊂ M i for all i.

Returning to the Michaelis-Menten example let us analyse the quasiequilibrium situation [START_REF] Meiske | An approximate solution of the michaelis-menten mechanism for quasi-steady and state quasi-equilibrium[END_REF][START_REF] Segel | On the validity of the steady state assumption of enzyme kinetics[END_REF][START_REF] Segel | The quasi-steady-state assumption: a case study in perturbation[END_REF][START_REF] Gorban | Asymptotology of chemical reaction networks[END_REF][START_REF] Gorban | The Michaelis-Menten-Stueckelberg theorem[END_REF] when the reaction constant k -1 is much faster than the reaction constant k 2 . In terms of orders, this condition reads γ -1 < γ 2 . In this case, the two tropical equilibration equations ( 9), ( 10

) are identical, because min(γ -1 , γ 2 ) = γ -1 . Let γ m = γ -1 -γ 1 denote the order of the parameter K m = k -1 /k 1 .
There are two branches of solutions of ( 9), namely a 2 = γ e , a 1 ≤ γ m and a 2 = a 1 + γ e -γ m , a 1 ≥ γ m corresponding to min(γ 1 + a 1 , γ -1 ) = γ 1 + a 1 and to min(γ 1 + a 1 , γ -1 ) = γ -1 , respectively. Using the relation between orders and concentrations we identify the first branch of solutions with the saturation regime x 2 ≈ e 0 (the free enzyme is negligible) and x 1 >> K m (the substrate has large concentration) and the second branch with the linear regime x 2 << e 0 (the concentration of the attached enzyme is negligible) and x 1 << K m (the substrate has low concentration).

The fast truncated system (obtained after removing all dominated monomials from ( 7)

) reads ẋ1 = -k 1 x 1 e 0 + k -1 x 2 , ẋ2 = k 1 x 1 e 0 -k -1 x 2 , ( 14
)
for the linear regime branch and ẋ1

= -k 1 x 1 (e 0 -x 2 ), ẋ2 = k 1 x 1 (e 0 -x 2 ), (15)
for the saturated regime branch.

Learning a finite state machine from a nonlinear biochemical network

Once the metastable regimes are computed, one can try to reconstruct the transitions happening between these regimes as illustrated in Fig. 1. By abstraction, we can consider that the metastable regimes (patches in the "crazy-quilt" picture in Fig. 1) are represented as nodes of a graph. Two nodes in the graph are connected if and only if there is at least one transition from one regime to the other. We call this abstraction a finite state machine, because the number of regimes is finite. However, given one regime, there may be several different possibilities to leave this regime, each leading to a different metastable regime. Which transition is chosen depends on the initial data. Thus, although the initial model is deterministic the finite state machine abstraction is generally stochastic. From one node can leave several transitions each one having a different probability per unit time. Because the initial system spends a long time on the metastable regime and little time on the transitions, it is natural to expect that the memory of previous transitions is lost and that the stochastic process is Markovian. In the following we propose a method to learn this stochastic, Markov process, from many simulations of the full model, each one starting from a different, randomly chosen, initial state. The method is similar to the Markov state models [START_REF] Bowman | Progress and challenges in the automated construction of markov state models for full protein systems[END_REF] used to coarse grain phase space for protein folding molecular dynamics. Contrary to energy landscape in protein folding, our "crazy-quilt" is not that rugged as it consists in a small number of patches. Therefore, standard Monte Carlo procedures function well in our case.

In order to construct the finite state machine, we first need a way to map the phase space of the continuous model to a finite set of branches. First, we compute the branches of tropical solutions as subsets of the euclidian space R n where n is the number of variables. We are using the algorithm based on constraint solving introduced by [START_REF] Soliman | A constraint solving approach to model reduction by tropical equilibration[END_REF] to obtain all rational tropical equilibration solutions a = (a 1 , a 2 , . . . , a n ) within a box |a i | < b, b > 0 and with denominators smaller than a fixed value d, a i = p i /q, p i , q are positive integers, q < d. The output of the algorithm is a matrix containing all the tropical equilibrations within the defined bounds. A post-processing treatment is applied to this output consisting in computing truncated systems, index sets, and minimal branches. Tropical equilibrations minimal branches are stored as matrices A 1 , A 2 , . . . , A b , whose lines are tropical solutions within the same branch. Here b is the number of minimal branches.

Our method computes numerical approximations of the tropical prevariety. Given a value of ε, this approximation is better when the denominator bound d is high. At fixed d, the dependence of the precision on ε follows more intricate rules dictated by Diophantine approximations. For this reason, we systematically test that the number b and the truncated systems corresponding to minimal branches are robust when changing the value of ε.

Trajectories x(t) = (x 1 (t), . . . , x n (t)) of the smooth dynamical system are generated with different initial conditions, chosen uniformly. For each time t, we compute the Euclidian distance d i (t) = min y∈A i y -log ε (x(t)) , where * denotes the Euclidean norm and log ε (x) = (log x 1 / log(ε), . . . , log x n / log(ε)). This distance classifies all points of the trajectory as belonging to a tropical minimal branch. The result is a symbolic trajectory s 1 , s 2 , . . . where the symbols s i belong to the set of minimal branches. In order to include the possibility of transition regions we include an unique symbol t to represent the situations when the minimal distance is larger than a fixed threshold. The choice of this threshold is robust (see discussion in Section 3.2 and Fig. 4). We also store the residence times τ 1 , τ 2 , . . . that represent the time spent in each of the state.

Special care should be taken when the model has a number of conservations laws. A conservation law is a linear combination of species concentrations that is kept constant during the dynamics, in other words an equation of the type

c i 1 x 1 + c i 2 x 2 + . . . + c i n x n = K i .
We suppose that the conservation laws are semi-positive (all c i are positive or nought). Then, several such equations together with the positivity conditions x i ≥ 0 define a polyhedron. We want to pick initial conditions for the trajectories x(t) uniformly in such a polyhedron. To start, each component j can take its initial concentration in the range [0, min(K i /c i j , such that c i j > 0)], but as we set the initial concentration of one component, the range available for components which are involved in the same constraint is reduced. The last component picked for each constraint must take the maximal value. Thus the sampling depends on the order in which initial concentrations are selected. To avoid introducing a bias related to this order, a random ordering of components is selected for each random initial state. Furthermore, if a component is picked as last in a constraint, its value is enforced (it must take all what remains). Conflicts may arise if this component is part of other constraints as well. To avoid this, we ensure that the last assigned item is specific to the constraint (not involved in any other constraint). This step may create problems with highly interdependent sets of constraints or constraints with less than two specific components. For the models considered here, this sampling works well as constraints have more specific components than overlapping ones.

By this method we generate N symbolic trajectories of length M defin-ing the vectors of successive states (s j 1 , s j 2 , . . . , s j M ) and residence times (τ j 1 , τ j 2 , . . . , τ j M ), where j ∈ [1, N ]. The stochastic automaton is learned as a homogenous, finite states, continuous time Markov process, defined by the lifetime (mean sojourn time) of each state T i , 1 ≤ i ≤ b and by the transition probabilities p i,j from a state i to another state j. We use the following estimators for the lifetimes and for the transition probabilities:

T i = ( N n=1 M m=1 τ n 1 s n m =i )/( N n=1 M m=1 1 s n m =i ) ( 16
)
p i,j = ( N n=1 M m=1 1 s n m =i,s n m+1 =j )/( N n=1 M m=1 1 s n m =i ), i = j, (17)
where 1 C stands for the indicator function equal to one if condition C is fulfilled or else equal to nought.

Results and discussion

3.1. Benchmarking metastable branch computation on Biomodels database Data source. For benchmarking, we selected 36 models from the r25 version of Biomodels database [START_REF] Le Novere | BioModels database: a free, centralized database of curated, published, quantitative kinetic models of biochemical and cellular systems[END_REF]. We have considered all the models that have polynomial vector field (mass action kinetics) and satisfy some technical constraints imposed by our Systems Biology Markup Language (SBML) parser (no function definitions, for instance). Models with zero valued parameters or which did not have at least one positive and one negative monomial per ODE were also filtered out (there were two of those).

Computation of minimal branches. The model SBML files are parsed and the polynomial vector fields are extracted. Thereafter, the conservation laws (that are the sum of the variables whose total concentration is invariant) are computed. The vector field along with the conservation laws are the input to an unpublished version of the tropical geometry based algorithm by Samal et al. (2015a) to compute the minimal branches. It should be noted here that due to the conservation laws the number of equations may exceed the number of chemical species.

According to the Eq.( 8) and to the geometric interpretation of tropical equilibrations from Sect. 2.1 the tropical solutions are either isolated points or bounded or unbounded polyhedra. Changing the parameter ε is just a way to approximate the position of these points and polyhedra by lattices or in other words by integer coefficients vectors. Finding the value of ε that provides the best approximation is a complicated problem in Diophantine approximation. For that reason, we preferred an experimental approach consisting in choosing several values of ε and checking the robustness of the results.

A summary of the analysis is presented in Table 1 with seven different choices for ε values (we wanted to have orders of magnitude close to decimal ones and to avoid commensurability between different values of ε; the choice 1/5, 1/7, 1/9, 1/11, 1/17, 1/19, 1/23 seemed good enough for this purpose). In addition, in Fig. 2 we plot the number of minimal branches versus the number of equations in the model. As can be noticed, this number is much lower than the number of states of a Boolean network with the same number of variables, which illustrates the advantage of our coarse graining with respect to other methods that discretize the values of the variables in order to obtain Boolean or multi-value networks [START_REF] Davidich | The transition from differential equations to boolean networks: a case study in simplifying a regulatory network model[END_REF].

The complete list of models and the corresponding statistics is provided as Supplementary File 1.

There are a few theoretical bounds for the number of tropical branches and metastable regimes. For monomolecular networks (linear differential equations) there can be at most one branch of full tropical equilibrations. However, if we allow that slow variables are not equilibrated in metastable regimes (cf. the notion of partial equilibration in Section 2.1) we showed elsewhere (Radulescu et al., 2015a) that there can be more than one metastable regime, but less than n, where n is the number of species. For nonlinear networks, one upper bound can be obtained by the tropical Bézout theorem [START_REF] Richter-Gebert | First steps in tropical geometry[END_REF], which says that the number of tropical equilibrations is exponential in n. However, our numerical experiments show that this bound is largely overestimated.

Table . 1 here Fig. 2 here

Biological significance of metastable branches for TGF-β signalling

As a case study we consider a nonlinear model of dynamic regulation of Transforming Growth Factor beta TGF-β signalling pathway that we have recently described by [START_REF] Andrieux | Dynamic regulation of Tgf-B signaling by Tif1γ: a computational approach[END_REF]. TGF-β signalling occurs through association of the ligand with TGF-beta type I (TGFBR1) and type II (TGFBR2) serine/threonine kinase receptors. TGF-β binding to TGFBR2 induces recruitment and phosphorylation of TGFBR1, which in turn transmits the signal through phosphorylation of SMAD2 transcription factor. Once phosphorylated, the SMAD2 hetero-dimerizes with SMAD4 and the complexes then migrate to the nucleus, where they regulate the transcription of TGF-β-target genes. In that context, the Transcriptional Intermediary Factor 1, TIF1-γ have been shown to function either as a transcriptional repressor or as an alternative transcription factor that promote TGF-β signalling. The apparent controversial effect of TIF1-γ on regulation of the SMAD-dependent TGF-β signalling was solved by a model integrating a ternary complex associating TIF1-γ with SMAD2 and SMAD4 complexes. This model has a dynamics defined by n = 18 polynomial differential equations and 25 biochemical reactions (the full set of ordinary differential equations can be found in the Appendix 1, the reaction scheme can be found in by [START_REF] Andrieux | Dynamic regulation of Tgf-B signaling by Tif1γ: a computational approach[END_REF] or in the SBML model that we provide as Supplementary File 2. The computation of the tropical equilibrations for this model shows that there are 9 minimal branches of full equilibrations (in these tropical solutions all variables are equilibrated). The connectivity graph of these branches (defined in Section 2.1) and the learned finite-state automaton (obtained with the method of Section 2.2) are shown in Fig. 3. The Table2 illustrates the convergence of the estimated transition probabilities when N , the number of Monte-Carlo samples, is increased. In order to illustrate the robustness of the classification of the states in the trajectories for each branch B i , i ∈ [1, 9] we computed the distribution of probability of Euclidian distances between randomly chosen states of the model, compatible with the conservation laws, and the branches. In Fig. 4 these distributions are compared with the distribution of minimal distances used to classify states on the model's trajectories. The latter distances are smaller and clearly separated from the random states distances.

Fig. 3 here Table 2 here Fig. 4 here The transition probabilities of the automaton are coarse grained properties of the statistical ensemble of trajectories for different initial conditions (cf. Section 2.1). Given a state and a minimal branch close to it, it will depend on the actual trajectory to which other branch the system will be close to next. However, when initial data and the full trajectory are not known, the automaton will provide estimates of where we go next and with which probability. For the example studied and for nominal parameter values, the branch B1 is a globally attractive sink: starting from anywhere, the automa-ton will reach B1 with probability one. This branch contains the unique stable steady state of the initial model. This calculation illustrates the basic properties of minimal branches of equilibrations. Trajectories of the dynamical system can be decomposed into segments that remain close to minimal branches. Furthermore, all the observed transitions between branches are contained in the connectivity graph resulting from the polyhedral complex of the tropical equilibration branches. This result proves the solidity of our tropical approach, because the geometric connectivity was not enforced to constrain the possible transitions; the fact that it is really the case emerges from the analysis of the trajectories. A change of parameter values can have several consequences: change the connectivity graph, change of the probabilities of transitions and change of the attractor position.

In order to understand the significance of the minimal branches and their relation with dynamic and physiologic properties of the network we have performed an analytic study of the tropical equilibration solutions. We show in the following (see also the Appendix 2) that the most important cause of the multiplicity of branches is the dynamics of the TGFBR1 and TGFBR2 receptors whose internalization and trafficking regulates TGF-β signalling (Le [START_REF] Roy | Clathrin-and non-clathrin-mediated endocytic regulation of cell signalling[END_REF]. These two receptors belong to a ligand-receptor module of 6 variables and 12 reactions that is decoupled from the rest of the network. More precisely, the ligand-receptor module activates the SMAD transcription factors but receives no feed-back (see Fig. 5) and can be studied independently from the rest of the variables. This module has been used with little variation in many models of TGF-β signalling [START_REF] Vilar | Signal processing in the TGF-beta superfamily ligand-receptor network[END_REF][START_REF] Zi | Constraint-based modeling and kinetic analysis of the Smad dependent TGF-beta signaling pathway[END_REF][START_REF] Chung | Quantitative modeling and analysis of the transforming growth factor β signaling pathway[END_REF].

Fig. 5 here We show in the Appendix 2 that the tropical equilibration of the ligand-receptor module form a two dimensional polyhedron conveniently parametrized by the concentration orders a 12 and a 13 of the receptors TGFBR1 and TGFBR2 respectively. The branches can be calculated analytically (Eqs.( 25) and ( 29) in Appendix 2). For the nominal values of the model parameters one of these branches is empty and the three remaining branches correspond to B 1 , B 2 and B 3 . The two other triplets of branches (B 4 , B 5 , B 6 ) and (B 7 , B 8 , B 9 ) correspond to the same mutual relations of variables in the ligand-receptor module. They are distinguished by the values of the remaining variables (the transcription factors module). Our computation of the automaton showed that the branches B i , i ∈ [4, 9] are practically inaccessible from states in branches B i , i ∈ [1, 3], therefore we will not discuss them here.

We have used symbolic computation to determine the steady states of the ligand-receptor module. This module has an unique steady state corresponding to concentrations orders that can be placed inside the polyhedron of tropical solutions using the Eq.( 4). The minimal branch containing the steady state is a sink of the coarse grained dynamics. The polyhedron of tropical solutions, its decomposition into minimal branches, and the position of the steady states inside it, depend on model parameters. Among model parameters two are important: k 18 and k 19 representing the production rate of the protein receptors TGFBR1 and TGFBR2, respectively. Consequently, these two parameters are correlated to gene expression and account for possible variability in mRNA levels of the two types of receptors. Fig. 6 shows the tropical equilibration branches of the ligand-receptor modules for various parameters k 19 corresponding to various TGFBR2 expression levels. For the nominal parameters used in the model, the branch B1 is a sink i.e. an attractor (the coarse-grained dynamics shows that the probability to leave this state is negligible), and the branches B2 and B3 are transient i.e. metastable. This means that starting in the branch B2 or B3 the receptor module will reach the branch B1 after a certain time and will remain there. However, over-expression of TGFBR2 modelled by changing the parameter k 19 and illustrated in the Fig. 6a-c can tilt the balance in favor of large concentrations of receptor of type 2 corresponding to branches B2, B3 of tropical solutions in the model. Interestingly, this change occurs by a displacement of the position of the steady state from B1 to B2 and B3 and not by a change of the concentration values allowed for these branches.

Fig. 6 here While [START_REF] Vilar | Signal processing in the TGF-beta superfamily ligand-receptor network[END_REF] have speculated that the ligand-receptor module is responsible for the versatility of the response of the TGF-β pathway, no experimental evidence support this hypothesis. Here, we now demonstrate that there are correlations between dynamical specificity characterized by membership to a particular branch of equilibration and cell phenotype. We illustrate such a comparison for the NCI-60 panel of cancer cell lines, a well established tool for tumor comparison and drug screening provided by the National Cancer Institute. Based on microarray analysis, these cell lines were found to cluster into two classes: epithelium-like (non-aggressive) and mesenchymal-like (aggressive) cell lines [START_REF] Ross | Systematic variation in gene expression patterns in human cancer cell lines[END_REF].

Using the global proteome analysis of this NCI-60 panel [START_REF] Gholami | Global proteome analysis of the nci-60 cell line panel[END_REF] we extracted the protein expression levels of TGFBR1 and TGFBR2 and showed that mesenchymal-like (aggressive) cell lines can be distinguished from epithelial-like (non-aggressive) cell lines by the increased level of TGFBR2 (Fig 6d.

The proteome data was compared to membership to tropical branches. According to Eq.( 4) there is a linear relation between opposite concentration orders -a i and logarithms of concentrations, -a i ≈ b log(x i ), i = 1, . . . , n (b = -1/ log(ε) > 0). We used opposite concentration orders -a i instead of a i because they change in the same direction as the concentrations (small opposite orders mean small concentrations and large opposite orders mean large concentrations). Therefore, in Fig 6a-c the relation T GF BR1 = T GF BR2 is verified on the bissector of the first quadrant, whereas T GF BR2 > T GF BR1 and T GF BR2 < T GF BR1 are valid above and below the bissector, respectively. When compared the proteome results with the membership to a particular branch of equilibration, we found that the distribution of concentration orders in branches place non-aggressive cancer cell lines in a range covered by branch 1, whereas the aggressive cancer cell lines are placed in a range covered by branches 2 and 3 (Fig 6d). Indeed, the ratio T GF BR2/T GF BR1 is small for the branch B 1 and in non-aggressive cancer cell lines, and is much larger for B 2 , B 3 and in aggressive cell lines.

Furthermore, we validated the association of up-regulation of TGFBR2 with mesenchymal-like appearance in an independent dataset of 51 breast cancer cell lines [START_REF] Neve | A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes[END_REF]. As we have recently described [START_REF] Ruff | The disintegrin and metalloprotease adam12 is associated with tgf-β-induced epithelial to mesenchymal transition[END_REF], comparative analyses between Basal B cell lines with mesenchymallike phenotype and Basal A and Luminal cell lines with epithelial morphology permitted to identify more than 600 differentially expressed genes that include TGFBR2. Gene expression data were now extracted for TGFBR1 and TGFBR2 and we showed that TGFBR2 gene expression is significantly induced in mesenchymal-like cell lines while TGFBR1 did not vary (Supplementary Fig. 1). In accordance with our observation, Parker et al. [START_REF] Parker | Lower expression levels of the transforming growth factor beta receptor type ii protein are associated with a less aggressive tumor phenotype and improved survival among patients with clear cell renal cell carcinoma[END_REF] have previously reported the association of low TGFBR2 expression with a lower aggressive tumor phenotype.

In summary, the tropical geometry analysis of the TGFβ signalling model first shows that the multiplicity of branches is due to the dynamics of TGFβ receptors. The more important parameters in this ligand-receptor module are the concentration of TGFBR1 and TGFBR2 and three main tropical branches are distinguished by the value of TGFBR2/TGFBR1 ratio (small in B1, intermediate in B2, large in B3). Importantly, we showed that the TGFBR2/TGFBR1 ratio is associated with tumor cell lines phenotype (high and low TGFBR2/TGFBR1 ratio for aggressive and non aggressive cell lines, respectively). Together these data demonstrated that tropical geometry analysis permits to discriminate between cellular states based on the evaluation of TGFβ receptors concentration. The importance of such up-regulation of TGFBR2 in aggressive cancer cell lines might be related to its implication in SMAD-independent signalling that includes PI3K-Akt, JNK, p38MAPK and Rho-like GTPases and which highly contribute to epithelial-mesenchymal transition [START_REF] Zhang | Non-Smad pathways in TGF-β signaling[END_REF][START_REF] Moustakas | Induction of epithelial-mesenchymal transition by transforming growth factor β[END_REF].

Together these observations suggest that metastable regimes defined by branches of minimal tropical equilibrations are associated with cell phenotypes. The idea of associating tropical minimal branches with clinical phenotype is similar to the idea of cancer attractors [START_REF] Huang | Cancer attractors: A systems view of tumors from a gene network dynamics and developmental perspective[END_REF] where the idea is that cancer cells are trapped in some abnormal attractors.

Conclusion

We have presented a method to coarse grain the dynamics of a smooth biochemical reaction network to a discrete symbolic dynamics of a finite state automaton. The coarse graining was obtained using a tropical geometry approach to compute the states. These states correspond to metastable dynamic regimes and to relatively slow segments of the system trajectories. The coarse grained model can be used for studying statistical properties of biochemical networks such as occurrence and stability of temporal patterns, recurrence, periodicity and attainability problems.

Further improvement and evolution is possible for this approach. First, the coarse graining can be performed in a hierarchical way. For the nonlinear example studied in the paper we computed only the full tropical equilibrations that stand for the lowest order in the hierarchy (coarsest model). As discussed in Section 2.1 we can also consider partial equilibrations when slow variables are not equilibrated and thus refine the automaton. Generally, there are more partial equilibrations than total equilibrations and learning an automaton on the augmented state set will produce refinements. Second, and most importantly, the dynamics within a branch could be also described. As shown elsewhere, reductions of the systems of ordinary differential equations are valid locally close to tropical equilibrations [START_REF] Noel | Tropical geometries and dynamics of biochemical networks application to hybrid cell cycle models[END_REF][START_REF] Noel | Tropicalization and tropical equilibration of chemical reactions[END_REF]Radulescu et al., 2015b;Samal et al., 2015b). Furthermore, the same reduction is valid for all the equilibrations in a branch. This suggests that a hybrid approach, combining reduced ODE dynamics within branch with discrete transitions between branches is feasible. The transitions can be autonomously and deterministically commanded by crossing the boundaries between branches that are perfectly determined by our approach.

The most important result of this paper is the extension of the notion of attractor to metastable regimes of chemical reaction networks and the proposition of a practical recipe to compute metastability. Metastable regimes correspond to low-dimensional hypersurfaces of the phase space, along which the dynamics is relatively slower. Most likely, metastable regimes have biological importance because the network spends most of its time in these states. The itinerancy of the network, described as the possibility of transitions from one metastable regime to another is paramount to the way neural networks compute, retrieve and use information [START_REF] Tsuda | Chaotic itinerancy as a dynamical basis of hermeneutics in brain and mind[END_REF] and can have similar role in biochemical networks. Our approach based on tropical geometry provides an algorithmic method to test these ideas further. The extension of this approach i.e. making use of statistical methods to compute the association of the tropical minimal branches with clinical phenotypes based on "-omics" data remains a topic of future research.

Property 4.2 (dominated first order reactions). If a species Y is consumed by several first order reactions of kinetic constants k 1 , k 2 , . . . , k r and if γ 1 ≤ γ 2 ≤ . . . ≤ γ p < γ p+1 ≤ γ p+2 ≤ . . . ≤ γ r , then the reactions k p+1 , . . . , k r can be eliminated and the resulting simplified tropical system has the same tropical equilibration solutions as the full system.

Proof of Property 4.2. The proof follows from the following obvious property of the min operation min(γ 1 , . . . , γ p , . . . , γ r ) = min(γ 1 , . . . , γ p ).

Using γ 26 < γ 20 , γ 28 < γ 21 (a condition satisfied by the nominal model parameters and meaning that internalization is more rapid than degradation for both receptors 1 and 2) and the Properties 4.1,4.2 we can justify the reduction illustrated in Fig. 7. Because the reduced model has the same tropical solutions as the full, larger model, it is enough to solve the tropical equilibration problem for the reduced model.

The polyhedron of tropical solutions defined by Eq.( 25) can be partitioned into minimal branches (also polyhedra). This can be done by checking which term is dominant in the ordinary differential equations for the variables x 12 , x 13 and x 14 (see Eqs. ( 18 1.0000 0.0000 0.9836 0.0164 0.0000 1.0000 0.0000 1.0000 0.0227 0.2727 0.0227 200 1.0000 0.0000 0.9675 0.0325 0.0000 1.0000 0.0000 1.0000 0.0217 0.2391 0.0109 500 1.0000 0.0000 0.9614 0.0386 0.0000 1.0000 0.0000 1.0000 0.0216 0.2381 0.0173 1000 1.0000 0.0000 0.9703 0.0264 0.0033 1.0000 0.0000 1.0000 0.0280 0.2258 0.0086 2000 1.0000 0.0000 0.9585 0.0365 0.0050 1.0000 0.0000 1.0000 0.0246 0.2173 0.0086 3000 1.0000 0.0000 0.9561 0.0389 0.0050 1.0000 0.0000 1.0000 0.0295 0.2342 0.0084 5000 0.9996 0.0004 0.9542 0.0399 0.0059 0.9996 0.0004 1.0000 0.0340 0.2460 0.0093 10000 0.9994 0.0006 0.9556 0.0363 0.0081 0.9973 0.0027 1.0000 0.0344 0.2370 0.0110

N 5→4 5→6 6→3 7→4 8→4 8→5 8→6 8→7 8→9 9→6 100
0.4091 0.2727 1.0000 1.0000 0.1373 0.4902 0.0784 0.0784 0.2157 1.0000 200 0.4130 0.3152 1.0000 1.0000 0.1300 0.5200 0.0600 0.0600 0.2300 1.0000 500 0.4156 0.3074 1.0000 1.0000 0.1098 0.5294 0.0784 0.0980 0.1843 1.0000 1000 0.4280 0.3097 1.0000 1.0000 0.1051 0.5333 0.0949 0.0909 0.1758 1.0000 2000 0.4497 0.2998 1.0000 1.0000 0.1121 0.5354 0.1020 0.0859 0.1646 1.0000 3000 0.4439 0.2840 1.0000 1.0000 0.1122 0.5315 0.0984 0.0951 0.1627 1.0000 5000 0.4227 0.2880 1.0000 1.0000 0.1160 0.5272 0.0969 0.0965 0.1633 1.0000 10000 0.4304 0.2872 1.0000 1.0000 0.1132 0.5344 0.0997 0.0932 0.1595 1.0000 Figure 1: Abstract representation of metastability as itinerant trajectory in a patchy phase space landscape. Dominant vector fields (red arrows) confine the trajectory to low dimensional patches on which act weak uncompensated vector fields (blue arrows). A typical trajectory contains slow segments within patches where dominant vector fields cancel, and transitions between patches in the fast direction of uncancelled dominant vector fields. Continuous (but non smooth) connections are also possible, corresponding to role reversal between dominant and dominated vector fields. The term crazy-quilt was coined to describe such a patchy landscape [START_REF] Gorban | Dynamic and static limitation in reaction networks, revisited[END_REF]. The different variables mean: x 12 : RI (TGBR1), x 13 : RII (TGFBR2), x 14 : LR (ligandreceptor complex), x 15 : RIe (TGFBR1 in endosome), x 16 : RIIe (TGFBR2 in endosome),

x 11 : LRe (LR in endosome).

Figure 2 :Figure 3 :Figure 4 :Figure 5 :Figure 6 :

 23456 Figure 2: Semi-log plot of the minimal branches versus the number of equations in the models from Biomodels repository for ε = 1/5. Comparison with a binary network number of states 2 n suggests sub-exponential scaling.

Figure 7 :

 7 Figure7: In order to compute the branches of tropical equilibration of the ligand-receptor module we use a reduced model. The reduced model is not necessarily a good approximation of the full dynamics, but has exactly the same tropical solutions as the full model. The different variables mean: x 12 : RI (TGBR1), x 13 : RII (TGFBR2), x 14 : LR (ligandreceptor complex), x 15 : RIe (TGFBR1 in endosome), x 16 : RIIe (TGFBR2 in endosome), x 11 : LRe (LR in endosome).

  This reads min(γ 18 , a 14 + γ 23 , a 12 + γ 26 ) = min(a 12 + a 13 + γ 22 + γ 35 , a 12 + γ 26 ) (19) min(γ 19 , a 14 + γ 23 , a 13 + γ 28 ) = min(a 12 + a 13 + γ 22 + γ 35 , γ 28 + a 13 ) (20) min(γ 24 , γ 25 , γ 23 ) + a 14 = a 12 + a 13 + γ 22 + γ 35 (21) Suppose now that the following condition is true min(γ 24 , γ 25 , γ 23 ) = γ 23 . (22) This condition is satisfied by the nominal parameters and, like the previous condition, means that receptors have relatively large life-times. Then from (21) we got a 14 = a 12 + a 13 + γ 22 + γ 35 -γ 23 and the equations (19),(20) become min(γ 18 , a 14 + γ 23 , a 12 + γ 26 ) = min(a 14 + γ 23 , a 12 + γ 26 ) (23) min(γ 19 , a 14 + γ 23 , a 13 + γ 28 ) = min(a 14 + γ 23 , a 13 + γ 28 ) (24) The solutions of (23), (24) can be easily found and form the following polyhedron ({a 12 + a 13 + γ 22 + γ 35 ≤ γ 18 } ∪ {γ 26 + a 12 ≤ γ 18 }) ∩ ({a 12 + a 13 + γ 22 + γ 35 ≤ γ 19 } ∪ {γ 28 + a 13 ≤ γ 19 }), a 14 = a 12 + a 13 + γ 22 + γ 35 -γ 23 . (25) The orders of the remaining variables can be found as indicated in Prop. 4.1: a 15 = a 12 + γ 26 -γ 27 , (26) a 16 = a 13 + γ 28 -γ 29 , (27) a 11 = a 12 + a 13 + γ 22 + γ 35 -γ 30 .

Table 2 :

 2 )). The result is that there are at most four minimal branches defined by one of the conditions {a 12 + γ 26 < a 12 + a 13 + γ 22 + γ 35 } ∩ {a 13 + γ 28 < a 12 + a 13 + γ 22 + γ 35 } {a 12 + γ 26 < a 12 + a 13 + γ 22 + γ 35 } ∩ {a 13 + γ 28 > a 12 + a 13 + γ 22 + γ 35 } {a 12 + γ 26 > a 12 + a 13 + γ 22 + γ 35 } ∩ {a 13 + γ 28 < a 12 + a 13 + γ 22 + γ 35 } {a 12 + γ 26 > a 12 + a 13 + γ 22 + γ 35 } ∩ {a 13 + γ 28 > a 12 + a 13 + γ 22 + γ 35 }(29) Fig. 7 here Estimation of the transition probabilities between branches of the TGFβ model, for different values of N , the number of Monte-Carlo samples.

	Table 1: Summary of analysis on Biomodels database. The benchmarked models have a
	number of dimensions (i.e. number of variables along with number of conservation laws)
	ranging from 2 to 41. Model BIOMD0000000289 has tropical branches at ε value 1/5, 1/7,
	1/9, 1/11 but none at 1/17, 1/19, 1/23.
	ε	Total	Models	Models	Average	Average	Min	Max
	value	mod-	without	with	running	num-	num-	num-
		els	tropical	tropical	time (in	ber of	ber of	ber of
		con-	branches	branches	secs)	minimal	minimal	minimal
		sidered				branches	branches	branches
	1/5	36	0	36	200.76	15.08	1	423
	1/7	36	0	36	177.01	14.41	1	406
	1/9	36	0	36	195.01	13.02	1	340
	1/11	36	0	36	169.59	12.36	1	322
	1/17	36	1	35	175.56	11.02	0	287
	1/19	36	1	35	187.42	11.08	0	287
	1/23	36	1	35	184.80	11.08	0	287

Acknowledgements. O.R and A.N are supported by INCa/Plan Cancer grant N • ASC14021FSA. D.G. is grateful to the grant RSF 16-11-10075 and to MCCME for wonderful working conditions and inspiring atmosphere.

Appendix 1: Description of the TGFb model used in this paper.

The model is described by the following system of differential equations dx

These variables are as follows:

• Receptors on plasma membrane: x 12 = RI (receptor 1), x 13 = RII (receptor 2), x 14 = LR (ligand-receptor complex).

• Receptors in the endosome: x 11 = LRe, x 15 = RIe, x 16 = RIIe.

• Transcription factors and complexes in cytosol:

• Transcription factors and complexes in the nucleus:

Appendix 2: Calculation of tropical equilibration branches for the ligand-receptor module.

Tropical equilibration solutions for the variables x 11 ,x 12 , x 13 , x 14 , x 15 , x 16 (the submodel in Fig. 5) can be computed independently from the rest of the variables of the TGFβ model. The ordinary differential equations for these variables form a subsystem that is decoupled (receives no feed-back) from the rest of the equations.

We can reduce the system of 6 tropical equations to a simplified system of 3 tropical equations using the following two general properties.

Property 4.1 (binomial species). Y is a binomial species if the ordinary differential equation defining its rate of variation contains only one positive monomial term and only one negative monomial term

where X denotes the other variables. We further assume that n 1 < n 2 . Then, the species Y can be eliminated and the resulting simplified tropical system has the same tropical equilibration solutions as the full system. The simplification is performed by eliminating the equation for Y and replacing everywhere Y by (M 1 /M 2 ) 1/(n 2 -n 1 ) .

Proof of Property 4.1. The proof follows from the fact that the tropical equation for the order a of Y has the unique solution a = 1 (n 2 -n 1 ) (µ 1 -µ 2 ).

Supplementary material captions.

Supplementary File 1. The complete list of models from the Biomodels database with the results of the tropical branches calculations.

Supplementary File 2. SBML version of the TGFβ model used in this study.

Supplementary Fig. 1. Gene expression levels for TGFBR2 and TGFBR1 in breast cancer cell lines from [START_REF] Neve | A collection of breast cancer cell lines for the study of functionally distinct cancer subtypes[END_REF]).


T Batal

F Ferlay

S Burles

S Larroque

L Doceul

F Samaille

M Missirlian

J Bucalossi

Mechanical Analysis of WEST divertor support plate

Keywords: WEST, Divertor, Plasma Facing Components, X-point configuration, actively-cooled

The Tore Supra tokamak is being transformed in an x-point divertor fusion device in the frame of the WEST (W-for tungsten-Environment in Steady-state Tokamak) project, launched in support to the ITER tungsten divertor strategy. The WEST project aims to test W monoblock Plasma Facing Units (PFU) under long plasma discharge (up to 1000s), with thermal loads of the same magnitude as those expected for ITER. Therefore the divertor is a key component of the WEST project, and so is its support structure, which has to handle strong mechanical loads. The WEST upper and lower divertor are made of 12 30° sectors, each one composed of 38 PFU that can be made of tungsten, CuCrZr or graphite. A generic 316L stainless steel 30° conic support plate is used to hold the 38 PFU together, regardless of their material. The PFUs are fixed on the support plate thanks to 152 Xm19 stainless steel fixing elements (4 per PFU), and in each of this fixing element an Aluminium-Nickel-Bronze alloy (Al-Ni-Br) pin is engaged in a slotted hole, in order to allow thermal expansion in the length direction of the PFU. The support plate is fixed on the divertor coil casing thanks to 10 M10 screws. Mechanicals loads which act on the PFUs are transmitted to the support plate through the fixing elements. These loads are due to Vertical Displacement Event (VDE), disruptions and thermal expansion of the PFU. First the different load cases, PFU configurations and scenario are presented. Then an ANSYS plastic mechanical simulation is performed in order to validate the number of cycles of the support plate for each scenario: 30 000 cycles in steady-state and 3000 cycles in VDE. Finally reactions forces from the previous ANSYS simulation are used in order to calculate the stress in the M10 screws.

Introduction

Tore supra is a tokamak using superconducting magnets. It was designed to carry out long plasma operations, so it uses actively cooled plasma facing components. The WEST project aims to transform Tore Supra from a limiter configuration to an X-point configuration in order to validate the technology of a full tungsten ITER-like actively cooled divertor [START_REF] Bucalossi | The WEST project: testing ITER divertor high heat flux component technology in a steady state tokamak environment[END_REF].

The W monoblock Plasma Facing Units (PFU) will be tested under long plasma discharge (up to 1000s), with thermal loads of the same magnitude as those expected for ITER. Therefore the divertor is a key component of the WEST project, and so is its support structure, which has to handle strong mechanical loads.

The WEST upper and lower divertor are made of 12 30° sectors, each one composed of 38 PFU that can be made of tungsten, CuCrZr or graphite. The graphite and tungsten PFU are assembled only on the lower divertor. The CuCrZr PFU can be assembled on lower and upper divertor. The divertor support plate has to be able to handle structural loads due to thermal expansion of the PFU and forces and torques resulting from halo and eddy current in the PFU during disruptions and Vertical Displacement Event (VDE, dBn/dt up to 90 T/s for 10 ms).

The WEST divertor support plate

The WEST divertor support plate is a generic 316L stainless steel 30° conic support plate. It is used indifferently for the lower and upper divertor, regardless of the material of the PFU. The PFUs are fixed on the support plate thanks to 152 Xm19 stainless steel fixing elements (4 per PFU), and in each of this fixing element an Aluminium-Nickel-Bronze alloy (Al-Ni-Br) pin is engaged in a slotted hole, in order to allow thermal expansion in the length direction of the PFU (Fig. 2). The support plate is fixed at its edges on the divertor coil casing thanks to 10 M10 screws. The screws are made of stainless steel BUMAX 109 (equivalent to class 10.9 bolts). In the inner part of the support plate, five of the M10 screws are engaged in the casing. In the outer part of the support plate, the other five M10 screws are engaged to supporting bracket welded to the casing. There is a 2 mm gap between the central part of the divertor support plate and the divertor coil casing (Fig. 1 and3).

Structural analysis

The structural analysis was performed with ANSYS 16.0.

Design criteria

The design criteria for the support plate is 30 000 cycles in steady-state and 3000 cycles in VDE. According to the ITER material handbook [START_REF]Appendix A, Materials Design Limit Data[END_REF], this corresponds to a maximal total strain (elastic strain+ plastic strain) for 316L stainless steel:  0,31% in steady states  0,54% in VDE The design criterion for the screws in steady-states is the bolt preload, which should not be overcome by external forces. The preload is 32 500 N for M10 class 10.9 bolt according to AFNOR standards (85% of yield stress).

As VDE is a very brief event, preload overcome can be allowed. The design criterion for the screws in VDE is the von Mises yield criterion which should not be higher than the yield stress, which is 900 MPa for BUMAX 109. The von Mises stress is defined by √𝜎 2 + 3𝜏 2 , where 𝜎 is the normal stress and 𝜏 the shear stress.

Boundary conditions

The temperature of the WEST support plate is 70°C. The mechanical behavior of the support plate is elastoplastic. The mechanical behavior of the screws is elastic, as the von Mises yield criterion has to be respected.

Surfaces under the head of the screws are fixed (Fig. 3). Compression supports are used to represent the divertor coils casing

Structural loads

As the fixing elements of the PFU allow displacement in their length direction, it is assumed that the forces acting on the PFU which are transmitted to the support plate are always locally perpendicular to its top surface. The structural loads acting on the divertor support plate are of two kinds:  Loads due to thermal expansion of PFU  Electromagnetic loads (EM loads) acting directly on the support plate or acting on the PFU and transmitted to the PFU through the fixing elements

The loads due to thermal expansion of the PFU are more significant with the upper divertor, which is composed of CuCrZr heatsinks with a tungsten layer. These thermal expansion loads were always considered in order to be conservative. These loads were studied in [START_REF] Batal | Thermomechanical simulation of WEST actively cooled upper divertor[END_REF]. These loads are the only ones in steady states.

During plasma disruption, the magnetic field inside the torus can vary rapidly. This will induce eddy current in the PFU or in the support plate. If the PFU assembled on the support plate are made of CuCrZr, the eddy current will first start to flow in the PFU because the electrical resistivity of CuCrZr (𝜌 = 2,10 • 10 -8 ohm.m at 70°C) is much lower than the one of 316L stainless steel (𝜌 = 7,40 • 10 -7 ohm.m at 70°C). So the PFU will generate a magnetic field opposite to the one that generated those eddy current, and will act as a magnetic screen for the support plate. If the PFU assembled on the support plate are made of graphite or tungsten monoblock, the eddy currents will first flow in the support plate which will act as a magnetic screen for the PFU. So during a plasma disruption, two configurations are possible:

If the PFU are made of CuCrZr, the eddy current flows in the PFU. The contribution of the magnetic field which is locally perpendicular to the PFU can be up to 90 T/s in upper position. This will induce a torque in the length direction of the PFU of 97 N.m [START_REF] Batal | Thermomechanical simulation of WEST actively cooled upper divertor[END_REF].

If the PFU are made of tungsten monoblock or graphite the eddy current flows in the support plate. The contribution of the magnetic field which is locally perpendicular to the support plate can be up to 20 T/s in lower position. By using an internal CEA ANSYS APDL routine [START_REF] Batal | Design and manufacturing of WEST Baffle[END_REF] an equivalent torque of 11 200 N.m was calculated (Fig. 6).

Fig 6. Disruption torque in the support plate

During a VDE the current flows through each PFU in the length direction because the plasma is in contact with the divertor [START_REF] Batal | Thermomechanical simulation of WEST actively cooled upper divertor[END_REF][START_REF] Larroque | The WEST project mechanical analysis of the divertor structure according to the nuclear construction code[END_REF]. Due to the toroidal contribution of magnetic field, which is 4.5 T, Halo forces are induced in the PFU, perpendicular to its top surface. This force is always directed from the top surface of PFU to its legs. According to [START_REF] Larroque | The WEST project mechanical analysis of the divertor structure according to the nuclear construction code[END_REF], the resulting force is -3,065 kN/PFU. Seismic loads have not been considered since TORE supra is not a nuclear installation. Dead weight and magnet fast discharge loads have not been considered either. So they are three different load cases in total that can be seen in Tab 1:

Steady

Structural Analysis

The simulations were run with ANSYS 16.0. The total strain in the support plate can be seen on Fig 7 to Steady-states 0,31 % 0,23 %  VDE 1 0,51 % 0,50 %  VDE 2 0,51 % 0,34 %  Tab. 2: Total strain in the three load cases

In order to check the stress in the screws, the reactions forces from the fixed support under the head of the screws are used. A normal force 𝐹 n and a radial force 𝐹 r are extracted from these reaction forces. A friction coefficient f of 0.5 is considered for the contact between the support plate on the divertor coil casings (in vacuum frictions forces are important). For the steady-states load case the equivalent preload 𝑃 𝑒 can be calculated thanks to the formula:

𝑃 𝑒 = 𝐹 n +
𝐹 r 𝑓 The equivalent section S of M10 screw is 58 mm². For the VDE load cases the von Mises stress 𝜎 𝑒 can be calculated thanks to the formula:

𝜎 𝑒 =√( 𝐹 n 𝑆 ) 2 + 3( 𝐹 r 𝑆 ) 2
The preload in steady-states and the von Mises stress in VDE are below the limits described in part 3.1. The results of the 10 screws can be seen in Tab. 3 and Tab. 4.

Screw

Conclusion

The structural integrity of the WEST divertor support plate was verified according the number of cycles defined in the WEST project: 30 000 cyclces in steadystates and 3000 cycles in VDE. The structural integrity of the screws fixing the support plate on the divertor coils casing was also verified: the preload was not overcome in steady-states and the von Mises yield criteria was respected in VDE. In this assessment the following load were considered: The most conservative thermal expansion loads of PFU expected during WEST operation in steady-states, and the eddy and halo current loads during disruptions and VDE. Seismic loads have not been considered since TORE supra is not a nuclear installation. Dead weight and magnet fast discharge loads have not been considered either.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The WEST upper divertor (radial-poloïdal plane)

  Fig 5.

Fig 4 .

 4 Fig 4. Fixed support under the heads of the screws

  Fig 9 and are summarized in Tab 2. The total strain in the three load cases is always below the limits defined in part 3.1.

Fig. 7 :

 7 Fig. 7: total strain in steady-states