E Gallet

M Manceny

P Le Gall

P Ballarini
email: paolo.ballarini@ecp.fr

An LTL Model Checking Approach for Biological Parameter Inference

Keywords: LTL Model Checking, Parameter Identification, Symbolic Execution, Genetic Regulatory Network, Thomas Discrete Modeling

The identification of biological parameters governing dynamics of Genetic Regulatory Networks (GRN) poses a problem of combinatorial explosion, since the possibilities of parameter instantiation are numerous even for small networks. In this paper, we propose to adapt LTL model checking algorithms to infer biological parameters from biological properties given as LTL formulas. In order to reduce the combinatorial explosion, we represent all the dynamics with one parametric model, so that all GRN dynamics simply result from all eligible parameter instantiations. LTL model checking algorithms are adapted by postponing the parameter instantiation as far as possible. Our approach is implemented within the SPuTNIk tool.

Introduction

Gene expression is a biological process where proteins are synthesized from genes. These proteins can regulate the synthesis of other proteins provided that their concentrations are sufficient. A collection of regulatory inter-dependencies between genes/proteins is called a Genetic Regulatory Network (GRN). In this paper we consider a discrete-state formalism, the René Thomas' formalism [START_REF] Barnat | On parameter synthesis by parallel model checking[END_REF][START_REF] Bernot | Application of formal methods to biological regulatory networks: extending Thomas' asynchronous logical approach with temporal logic[END_REF][START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Thieffry | Formalisation of regulatory networks : a logical method and its automation[END_REF][START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. immunity control in bacteriophage lambda[END_REF], according to which the amounts of proteins in a GRN are discrete abstractions of continuous concentrations. Hence the overall evolution of protein concentration along time, called the dynamics of the network, is captured by a discrete-state transition system. Given a René Thomas GRN model, the main interest is in analyzing the possible dynamics that may be associated to it. However, from the dynamics point of view, a GRN on its own is an underspecified type of representation: it represents the dependencies between a set of genes, but it does not describe what effect the combination of all such dependencies has on a given state of the network, hence on its evolution. The mapping of a GRN model to a specific dynamic (i.e. a transition system) is achieved by considering an instantiation of so-called biological parameters, i.e., a specification of the combined effect that all activated regulators have on a given state of the network.

Since the number of possible instantiations of biological parameters is a double exponential function of the GRN size (in terms of number of genes and of interactions), the analysis of the possible dynamics associated to a GRN model is a complex task. In this paper we consider the application of formal methods, namely model checking [START_REF] Baier | Principles of Model Checking[END_REF], as a means to reason about the possible dynamics associated to a GRN specification. In particular we tackle the following problem: given a relevant behavior of interest, formally expressed in terms of a temporal logic property, say ϕ, we want to be able to automatically identify the biological parameters instances which give rise to dynamics complying with ϕ.

Related work. Model checking techniques have been widely advocated in several works to verify whether a given discrete Thomas model fulfills some relevant biological temporal properties. In [START_REF] Bernot | Application of formal methods to biological regulatory networks: extending Thomas' asynchronous logical approach with temporal logic[END_REF] Bernot et al. expressed biological knowledge with Computation Tree Logic (CTL) formulas [START_REF] Baier | Principles of Model Checking[END_REF]. To exhaustively search the parameters' space, the set of all possible dynamics is generated and a CTL model checking procedure is iterated, one dynamics after the other. This approach is implemented in the SMBioNet tool [START_REF] Richard | SMBioNet User manual[END_REF] and has been illustrated in [START_REF] Filopon | Epigenetic acquisition of inducibility of type III cytotoxicity in P. aeruginosa[END_REF]. In [START_REF] Khalis | The SMBioNet method for discovering models of gene regulatory networks[END_REF], the approach has been extended to cope with the formation of complexes from proteins which allows modelers to express relationships between biological parameters leading to a reduced set of dynamics to be investigated. This work prefigures the interest of using constraints on the parameters. [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF][START_REF] Barnat | On parameter synthesis by parallel model checking[END_REF] define an approach based on an encoding technique to share computations between different dynamics. Sets of dynamics are encoded by a binary vector, one bit (or color) per dynamics, and LTL model checking algorithms [START_REF] Baier | Principles of Model Checking[END_REF] are extended with Boolean operations on vectors. In [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Corblin | Applications of a formal approach to decipher discrete genetic networks[END_REF], the tool GNBox uses Constraint Logic Programming (CLP) techniques to identify parameters. GRN dynamics and biological knowledge are described by declarative rules and constraints on parameters, then target behaviors are expressed as some kind of finite paths that models have to verify. [START_REF] Fromentin | Analysing gene regulatory networks by both constraint programming and model-checking[END_REF] also uses CLP techniques to adapt CTL model checking, but the encoding introduces a lot of fresh logical variables that hamper to scale up the method.

Our contribution. We propose a new approach that is based on a parametric model, called Parametric GRN (PGRN). This allows us to encompass all the dynamics of a GRN in a unique representation, biological parameters being processed as symbols, and to implement an efficient (on-the-fly) searching of (a symbolic representation of) the parameter's space. Similarly to [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF], our approach is based on LTL model checking. While in, [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF], LTL model checking algorithms are optimized for the particular case of time series, that are sequences of states made of one expression level per gene, observed one by one, we follow the same creed as the one advocated in [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Corblin | Applications of a formal approach to decipher discrete genetic networks[END_REF]: model sets are handled through some logical language both to avoid combinatorial explosion and to take benefit of constraint solving techniques. A preliminary version of our approach has been described in [START_REF] Mateus | Symbolic modeling of genetic regulatory networks[END_REF]. In the present version, algorithms combining symbolic execution and constraint solving techniques have been reengineered and tuned to be more efficient and cope with GRN features. Thereby, we consider the full LTL language while [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF] essentially focuses on time series and [START_REF] Corblin | A declarative constraint-based method for analyzing discrete genetic regulatory networks[END_REF][START_REF] Corblin | Applications of a formal approach to decipher discrete genetic networks[END_REF] focuses on properties carrying on finite paths.

Paper organization. We reformulate the logical description of Thomas' modeling framework in section 2, and explain how we encode the set of dynamics of a GRN with Parametric GRN in section 3. Section 4 presents our adaptation of LTL model checking algorithms with symbolic execution techniques. In section 5, we briefly discuss the validity of our approach with our dedicated tool SPuTNIk. Finally, section 6 contains some concluding remarks.

Genetic Regulatory Networks

A Genetic Regulatory Network (GRN) is a collection of regulatory inter-dependencies between genes to represent the mechanism of gene expression, i.e. the biological process by means of which proteins are synthesized. Two kinds of interactions exist: activation or inhibition depending on whether the protein expressed from the source gene can enhance or reduce the expression of the target gene. Moreover, an interaction is effective only if the concentration of the source protein is sufficient, in other words if its level of expression is above a given threshold. From the modeling point of view, a gene g is assimilated to the protein it synthesizes. In particular, it inherits the protein's level of expression, denoted x g , which, in this context, is abstractly represented by a non-negative integer ranging from 0 (absence of protein or very low concentration level) to a maximal value m g . A GRN is classically represented by an interaction graph.

Definition 1 (Interaction graph

). An n-order interaction graph (IG) is a labeled directed graph Γ = (G, I) where G is a finite set of gene nodes, n = |G| and I ⊆ G × {+, -} × N + × G is the set of interactions. Given (g 1 , s, t, g 2 ) ∈ I, s indicates the effect of g 1 over g 2 (sign "+" for activation and "-" for inhibition) and t denotes the threshold of the interaction. Moreover, the following properties hold: i) ∀(g 1 , g 2 ) ∈ G 2 , there exists at most one interaction (g 1 , s, t, g 2 ) ∈ I; ii) ∀(g 1 , s, t, g 2 ) ∈ I, t > 1 ⇒ ∃(g 1 , s , t , g 3 ) ∈ I, t = t -1.

The threshold t of an interaction (g 1 , s, t, g 2 ) ∈ I indicates the minimal level that g 1 needs to be at in order to affect the expression of g 2 . The condition on thresholds states that for any gene g 1 every intermediate threshold level must appear on at least an interaction arc originating in g 1 . Since an interaction graph may contain at most one interaction between two genes, then for an interaction (g 1 , s, t, g 2 ) ∈ I, we denote s(g 1 , g 2 ) its sign, and t(g 1 , g 2 ) its threshold. m g = max{t | ∃(g, s, t, g ) ∈ I} is 3 the maximal level of expression of gene g, and G -(g) ⊆ G is the set of regulators of g (i.e. G -(g) = {g |∃(g , s, t, g) ∈ I}).

A dynamics of a GRN corresponds to an evolution over time of the levels of expression of all genes, The state space describes the states that may be observed during such a possible evolution.

Definition 2 (State space of an interaction graph). For Γ = (G, I) an interaction graph we define X = g∈G X g the state space underlying Γ , where X g = {0, ..., m g } is the set of possible levels of expression for gene g.

Example 1 (Interaction graph and state space). Figure 1 presents a two-genes interaction graph Γ 0 where gene α is both an activator of β and of itself (selfactivator) while β is an inhibitor of α. In particular α activates the expression of β whenever its level of expression is at least 1, while when its level of expression is at least 2 it activates both itself and β. The thresholds of Γ 0 induce the following sets of levels for the two genes, X α = {0, 1, 2} and X β = {0, 1}, hence the state space X = {(0, 0), (0, 1), (1, 0), (1, 1), (2, 0), (2, 1)}.

α β [+, 2] [+, 1] [-, 1]
Fig. 1: Γ 0 : an example of a two-genes interaction graph.

xα x β (0, 0) (1, 0) (2, 0) (0, 1) (1, 1) (2, 1)
Fig. 2: Dynamics D K induced by the parameters mapping K

For Γ an interaction graph with state space X, we denote D = (X, →) a generic transition system4 called dymanics of Γ . With respect to D, Γ can be regarded as an underspecified formalism: from Γ , on its own, one cannot devise any concrete dynamics D. To obtain a dynamics for Γ we need to describe for each gene g the effect that any subset of its regulators ω ⊆ G -(g) would have on g. This is achieved by associating biological parameters with Γ . Definition 3 (Biological parameters). For g a gene of Γ = (G, I), the set of biological parameters of g is

K g = {K g (ω) | ω ⊆ 2 G -(g) } and the set of all biological parameters of Γ is K = ∪ g∈G K g .
An instantiation of biological parameters is defined by any mapping K : K → ∪ g∈G X g associating to any parameter K g (ω) a value in X g . Any instantiation K : K → ∪ g∈G X g defines a mapping K :

⊗ g∈G 2 G -(g) → X verifying ∀(ω g1 , . . . , ω gn ) ∈ ⊗ g∈G 2 G -(g) , K((ω g1 , . . . , ω gn )) = (K(K g1 (ω g1 )), . . . , K(K gn (ω gn ))) (with G = {g 1 , . . . , g n }).
In the sequel, for simplicity purpose, an instantiation K associating the value x to the parameter K g (ω) will be simply given by the equality K g (ω) = x. The biological parameters of an IG indicate the values the genes of the GRN tend towards when a certain n-tuple of regulators is activated. Thus to obtain the dynamics D K corresponding to parameters K we need to know what regulators are activated in a state x ∈ X. We say that a set of regulators ω ⊆ G -(g) of gene g is activated in state

x = (x 1 , . . . , x n ) ∈ X, denoted (x 1 , . . . , x n ) |= ω, iff g ∈ω x g ≥ t(g , g) ∧ g ∈G -(g)\ω x g < t(g , g)
, that is: if and only for every regulator g ∈ ω the corresponding component x g of state x is above the corresponding activation threshold (i.e., x g ≥ t(g , g)) while no other regulator of g does. In the remainder we denote ActR : X → ⊗ g∈G 2 G -(g) the function that maps each state x into the corresponding n-tuples of activated regulators for the n genes of a GRN. In the remainder, for x ∈ X, we denote x[x g ↑] (resp. x[x g ↓]) the state resulting from x by increasing (resp. decreasing) the x g component of one unit.

Definition 4 (dynamics induced by an instantiation of parameters).

For Γ = (G, I) an n-order IG, K : K → ∪ g∈G X g a set of biological parameters of Γ , we define D K = (X, → K ) the dynamics (transition system) of Γ induced by K. The transition relation → K ⊆ X ×X is minimally defined as follows:

∀x = (x 1 , . . . x n ) ∈ X let x * = (x * 1 , . . . , x * n ) = K • ActR((x 1 , . . . , x n )): -if x = x * then ∀i, 1 ≤ i ≤ n • if x i < x * i , then x → K x[x i ↑] (increment gene i) • if x i > x * i , then x → K x[x i ↓] (decrement gene i) -else if x = x * then x → K x (self-loop)
Then for each state x ∈ X we determine the corresponding attractor state x * (by application of the parameters mapping K to the regulators activated in x i.e., x * = K • ActR(x)). If the attractor state x * is different from x then for each different component x i = x * i we add either an increment or a decrement transition in → K . Conversely if x * = x we add a self-loop in → K .

Example 2 (Biological parameters and Dynamics). The subsets of regulators for the two genes of Γ 0 (Figure 1) are 2 G -(α) = {{}, {α}, {β}, {α, β}}, resp. 2 G -(β) = {{}, {α}} hence the biological parameters of Γ 0 are:

K = {K α ({}), K α ({α}), K α ({β}), K α ({α, β}), K β ({}), K β ({α})}.
According to |= the association between states of X and n-tuples of activated regulators is: (0, 0) |= ({}, {}), (0, 1) |= ({β}, {}), (1, 0) |= ({}, {α}), (1, 1) |= ({β}, {α}), (2, 0) |= ({α}, {α}), (2, 1) |= ({α, β}, {α}). As an example we consider the following mapping (instantiation) of the parameters for gene α and β into corresponding target levels: 2 shows the dynamics D K yielded by the parameters mapping K.

K α ({}) = 2, K α ({α}) = 2, K α ({β}) = 0, K α ({α, β}) = 2, K β ({}) = 0 and K β ({α}) = 1 yielding the combined mapping K = K α ×K β = {({}, {}) → (2, 0), ({{}, {α}) → (2, 1), . . . , ({α, β}, {α}) → (2, 1)} Figure
A parameters mapping K for an IG Γ yields a dynamics D K . However some mappings K may result into inconsistent dynamics. To rule out inconsistent dynamics, mapping must comply with the following constraints.

Definition 5 (Constraints for parameters mapping). Let

Γ = (G, I) be an IG. Definition constraint: ∀g ∈ G, ∀g ∈ G -(g), ∀ω ⊆ G -(g)\{g }: if s(g , g) = + then K g (ω) ≤ K g (ω∪{g }), if s(g , g) = -then K g (ω) ≥ K g (ω∪{g }). Observation constraint: ∀g ∈ G, ∀g ∈ G -(g), there exists ω ⊆ G -(g) \ {g }: if S(g , g) = + then K g (ω) < K g (ω ∪{g }), if s(g , g) = -then K g (ω) > K g (ω ∪{g }). Min/Max constraint: ∀g ∈ G, K g ({g |g ∈ G -(g), s(g , g) = -}) = 0 and K g ({g |g ∈ G -(g), s(g , g) = +}) = m g .
The Definition constraint (or Snoussi constraint [START_REF] Snoussi | Logical identification of all steady states: the concept of feedback loop characteristic states[END_REF]) states that if the level of expression of a gene g which activates (resp. inhibits) a gene g becomes greater than its threshold, then the expression level of g cannot decrease (resp. increase). The Observation constraint expresses how we identify regulators. If g is an activator (resp. inhibitor) of g, then there exists at least one dynamic state where the increase of the level of expression of g leads to an increase (resp. decrease) of the expression level of g. Finally, the Min/Max constraint states that in a dynamic state where all the activators (resp. inhibitors) of a gene are above the threshold and simultaneously none of the inhibitors (resp. activators) is, then the level of expression of the attractor of the gene is maximum (resp. minimum).

Example 3. The Constraints for IG Γ 0 (Figure 1) correspond to the following conditions:

K α ({α}) = 2, K α ({β}) = 0, K β ({}) = 0, K β ({α}) = 1, K α ({}) < 2 ∨ 0 < K α (α, β) and K α ({}) > 0 ∨ 2 > K α (α, β)
. Notice that amongst the 324 possible parameter mappings5 for Γ 0 , only 7 are consistent with the Constraints for Γ 0 .

Even if these constraints are well-founded, there are not always considered by biologists. In the sequel, by default, they will be considered and generically denoted as C I , but they can be relaxed on demand.

Modeling dynamics with Parametric GRN

In order to study all the dynamics simultaneously, we represent them all through a single (meta)model, called Parametric GRN (PGRN), i.e. a facility of transition systems parameterized by the biological parameters. Parametric GRN. A PGRN is a transition system associated with an interaction graph Γ = (G, I). It involves two families of symbols: the biological parameters

K = {K g (ω) | g ∈ G, ω ⊆ 2 G -(g)
} and the state variables G = {x g |g ∈ G}. Note that, according to the context, x g will denote either a state variable or a value representing a concentration level.

The main idea is to encode state evolution with transitions parameterized by parameters of K. A PGRN is composed of two states: T (transient) corresponding to configurations such that at least one gene can change its current level, and S (stable) corresponding to situations where no change is possible for any gene. A transition of a PGRN is characterized by a guard (a condition over parameters K and state variables of G) and an assignment (an application X → X expressing how states of a dynamics evolve). For example, transition T

(xα<2∧x β =0∧xα<Kα({})[xα↑]
-------------------→ T (see Fig. 3) indicates that for any (transient) state

x ∈ X such that x β = 0, x α < 2 and x α < K α ({}) then a transition corresponding to an increase of the level of α exists.

More precisely, for each gene g, there is a transition from T to T for each kind of variation (increase or decrease) of x g . For ω ⊆ G -(g) a subset of regulators of g, let us introduce the predicate P g (ω) : X → { , ⊥} defined by: ( g ∈ω x g ≥ t(g , g)) ∧ ( g ∈G -(g)\ω x g < t(g , g)). P g (ω) characterises the set of states in which regulators ω are the only effective ones on g. The transition associated to the increase of x g is conditioned by the guard Increase(g) = ∨ ω⊂G -(g) (P g (ω) ∧ x g < K g (ω)). Similarly, the transition associated to the decrease of x g is conditioned by the guard Decrease(g) = ∨ ω⊂G -(g) (P g (ω) ∧ x g > K g (ω)). Finally there is one transition from T to S when the expression level of all genes remains stable, i.e. if any gene g satisfies the condition Stable(g) = ∧ ω⊂G -(g) (P g (ω) ∧ x g = K g (ω)), and one last transition from S to S where the guard is always true.

Definition 6 (PGRN).

A PGRN associated to an interaction graph Γ = (G, I) is a pair P = (Q P , δ P ) with Q P = {T, S} the set of states and δ P a set of transitions. A transition of δ P is of the form (q P , g P , a P , q P ), also denoted

q P (g P )[a P ]
-----→ q P , with q P and q P states of Q P , g P a guard, i.e. a formula over K∪G and a P an assignment, i.e. an application X → X. More precisely, δ P is the set of all following transitions:

-(T, Increase(g), x g ↑, T ) with g in G, -(T, Decrease(g), x g ↓, T ) with g in G, -(T, ∧ g∈G Stable(g), id, S) where id is the identity assignment, -(S, , id, S) where indicates the guard always true.

Let us remark that unfolded versions of guards can be rather long and complex, but in the best cases they can be simplified by application of the initial constraints C I . Nevertheless, generally, the most complex guard is the one labeling the transition T → S since it corresponds to the conjunction of all Stable(g) conditions. On the other hand, once in S, the guard of the only possible transition (S → S) is simply true ( ). Moreover, transitions involving disjunctions in their guard can be split. Indeed, transition (T, g P ∨ g P , a P , T ) can be equivalently split in (T, g P , a P , T ) and (T, g P , a P , T ).

Example 4 (PGRN). Fig. 3 represents the PGRN associated with the interaction graph of Fig. 1. In relation with the different possible subsets ω, one can explicit the different guards: e.g.

Increase(β) ≡ (x α < 1∧x β < K β ({}))∨(x α ≥ 1∧x β < K β ({α}))
). The Initial constraints C I (cf. Def 5) can be used to simplify the guards: e.g. C I implies K β ({}) = 0 and K β ({}) = 1 and then Increase(β) ≡

(x α > 0 ∧ x β = 0). T S (xα < 2 ∧ x β = 0 ∧ xα < Kα({}))[xα ↑] (xα > 0 ∧ x β = 0)[x β ↑] (xα = 0 ∧ x β = 1)[x β ↓] (xα = 1 ∧ x β = 0 ∧ Kα({}) = 0) ∨ (xα = 1 ∧ x β = 1)∨ (xα = 2 ∧ x β = 1 ∧ Kα({α, β}) < 2)[xα ↓] (xα = 0 ∧ x β = 0 ∧ Kα({}) = 0)∨ (xα = 2 ∧ x β = 1 ∧ Kα({α, β}) = 2)[id]
[id]

Fig. 3: PGRN associated to the interaction graph in Figure 1.

Annotated dynamics. A PGRN characterizes a set of dynamics, one for each possible instantiation of biological parameters, that is, for any parameter mapping K : K → ∪ g∈G X g . For g P a transition guard of a PGRN, x ∈ X a state of the corresponding GRN, and K an instance of biological parameters, we denotes g P x,K the instance of g P obtained by substituting g P 's state variables and g P 's biological parameters with the corresponding state value of x, and parameter values of K. Similarly, we denote g P x the resulting substitution only of g P 's state variables (parameters in K remain symbolic).

Definition 7 (Annotated Dynamics). Let P = (Q P , δ P ) be a PGRN associated with an interaction graph Γ = (G, I), and let K : K → ∪ g∈G X g . The annotated dynamics associated to P and K is a pair

D K = (Q D , δ D )
where the set of states Q D ⊂ Q P × X and the set of transitions δ D ⊂ Q D × Q D are mutually defined by: ∀x ∈ X, (T, x) ∈ Q D and for all (q P , x) ∈ Q D and (q P , g P , a P , q P ) ∈ δ P s.t. g P x,K is evaluated to True, then (q P , a P (x)) ∈ Q D and ((q P , x), (q P , a P (x))) ∈ δ D .

Example 5 (Annotated Dynamics). Figure 4 presents one possible annotated dynamics for the PGRN represented in Figure 3, with the following instantiation of parameters:

K α ({}) = 2, K α ({α}) = 2, K α ({β}) = 0, K α ({α, β}) = 2, K β ({}) = 0 and K β ({α}) = 1. T, (0, 0) T, (1, 0) T, (2, 0)
T, (0, 1) T, (1, 1) T, (2, 1) S, (2, 1) Fig. 4: A possible annotated dynamics for the PGRN in Figure 3.

By construction, for a given instantiation K : K → ∪ g∈G X g , the associated annotated dynamics D K corresponds to the dynamics D K of the underlying IG Γ induced by the instantiation K (cf Def 4). For a transition ((q P , x), (q P , a P (x))) in δ D , it suffices to give up the first component and keep the second one, x → a P (x), to retrieve a dynamics of Γ . Thus, the dynamics represented in Fig. 2 can be obtained from the annotated dynamics of Fig. 4. The first component (T or S) is somehow a technical artifact annotating the presence of a stable state when building sequences of consecutive states. Depending on the context, we will assimilate D K and D K or will work with the most appropriate of the two forms. Motivated by efficiency considerations, we will apply a specific treatment for states x already recognized as stable, that is, annotated by S.

Adapting LTL model-checking to PGRN

The classical approach of LTL model-checking [START_REF] Pnueli | The temporal logic of programs[END_REF] consists in confronting a model (e.g. a dynamics) against an LTL formula. To do so, the negation of the LTL formula is transformed into a Büchi automaton and the product between the automaton and the dynamics is computed. We then look for accepting paths in the product by checking the existence of reachable cycles containing at least an accepting state. Model checking is usually time consuming, and since the number of dynamics is large, this method is not applicable in our case. To avoid the combinatorial explosion, we want to check all the dynamics simultaneously, i.e. we check directly the PGRN. To do so, we first build the Parametric Product between the PGRN and the Büchi Automaton associated to the LTL formula ϕ. We then use symbolic execution technics in order to search for accepting cycles. As a result, we obtain a set of constraints C that a parameter instantiation K must fulfill such that the associated dynamics D K verifies ϕ.

Büchi Automaton and Parametric product Biological properties on a sequence of states can be expressed using LTL formulas built from a set of atomic propositions using the usual logical operators in { , ⊥, ¬, ∧, ∨} and the temporal operators X (for neXt time), G (Globally), F (Finally) and U (Until) [START_REF] Baier | Principles of Model Checking[END_REF]. Since we need to express biological knowledge on levels of expression of genes, atomic propositions are of the form x g c where x g denotes the level of expression of a gene g, ∈ {=, =, <, >, ≤, ≥} and c ∈ N. Any LTL formula ϕ can be translated into a Büchi automaton B(ϕ). Definition 8 (Büchi Automaton associated to an LTL formula). Let Γ be a GRN and ϕ an LTL formula over the levels of expression of genes of Γ . A Büchi Automaton associated to ϕ is a tuple

B(ϕ) = (Q B , q 0 B , A B , δ B ) where Q B is the set of states, q 0 B ∈ Q B is the initial state, A B ⊆ Q B is
the set of accepting states and δ B is the set of transitions. A transition of δ B is of the form (q B , g B , q B ) with q B and q B states of Q B and g B a non temporal formula over the levels of expressions of genes in Γ . Moreover, B(ϕ) is such that an infinite sequence of states provided with truth values for all atomic propositions (a path) verifies ϕ iff this path is accepted by B(ϕ), i.e. iff this path contains at least a so-called accepting state infinitely often.

Example 6 (LTL formula and associated Büchi automaton). The existence of a steady state (i.e. a state which is itself its only own successor) in (x α , x β ) = (2, 1) corresponds to the LTL formula G((x α = 2 ∧ x β = 1) → X(x α = 2 ∧ x β = 1)). Fig. 5 presents a Büchi Automaton associated to the negation of this formula.

s0 start s1 s2 s3 xα = 2 ∧ x β = 1 xα = 2 ∧ x β = 1 xα < 2 x β = 0 Fig. 5: B(¬ϕ) with ϕ ≡ G (x α = 2 ∧ x β = 1) ⇒ X(x α = 2 ∧ x β = 1) .
Definition 9 (Parametric Product). Let P = (Q P , δ P ) be a PGRN and

B(¬ϕ) = (Q B , q 0 B , A B , δ B ) a Büchi Automaton associated to the LTL formula ¬ϕ. The product Π = P ⊗ B(¬ϕ) is the tuple (Q Π , q 0 Π , A Π , δ Π ) with Q Π = Q P × Q B the set of vertices, q 0 Π = (T, q 0 B ) the initial vertex, A Π = Q P × A B the
set of accepting vertices, and δ Π the set of transitions. A transition of δ Π is of the form (q Π , g Π , a P , q Π ) with q Π = (q P , q B ), q Π = (q P , q B ), g Π = g P ∧ g B such that (q P , g P , a P , q P ) ∈ δ P , (q B , g B , q B ) ∈ δ B and g Π is satisfiable.

Example 7. The product of the PGRN in Fig. 3 and the Büchi Automaton in Fig. 5 is represented in Fig. 6. The product has been simplified by removing output transitions whose guard on expression levels is not satisfiable according to the guards and assignments of the input transitions of the same vertex; we also remove the transitions whose guard is not satisfiable according to the guards on parameters necessarily crossed (φ 21 and φ 22 here). Finally, we remove vertices which can not be reached and those belonging to a terminal cycle without accepting vertex.

Search for parametric accepting cycles The search for accepting cycles is based on symbolic execution techniques which are program analysis techniques.

The key point is the substitution of actual values by symbolic variables in order to symbolically perform computations. Each execution (or path) of the program associates to each variable a symbolic computation together with a path condition that expresses what are the conditions on input values to execute the given path. Symbolic execution techniques has been extended to symbolic transition systems [START_REF] Gaston | Symbolic execution techniques for test purpose definition[END_REF] by unfolding transition systems as symbolic trees. As symbolic execution is only applicable for finite paths, selection criteria are used to cut infinite paths when considering testing. In the sequel, we will take particular care to cut infinite paths in identifying situations of return on a node already encountered. Indeed such situations reveal the presence of cycles.

In the symbolic execution of the parametric product Π, the parameters K g (ω) are handled as symbolic variables (i.e. not evaluated), and Π is unfolded leading to the construction of several Symbolic Execution Trees (SET), one for any x ∈ X.

Definition 10 (Symbolic Execution Tree). Let Π = (Q Π , q 0 Π , A Π , δ Π ) be a parametric product. The Symbolic Execution Tree associated to Π and x ∈ X is a transition system (Q T , δ T ) where the set of nodes Q T and the set of transitions

Ts0 start Ts1 Ts3 Ss3 φ1[xα ↑] φ2 1 [xα ↓] φ2 2 [xα ↓] φ2 3 [xα ↓] φ3[x β ↓] φ4[x β ↑] φ2 1 [xα ↓] φ2 2 [xα ↓] φ1[xα ↑] φ2 1 [xα ↓] φ2 2 [xα ↓] φ2 3 [xα ↓] φ3[x β ↓] φ4[x β ↑] φ5[id]
[id]

Fig. 6: Parametric product P ⊗ B(¬ϕ) associated to the Parametric GRN in Figure 3 and the Büchi Automaton in Figure 5 (after simplification), with:

φ1 ≡ xα < 2 ∧ x β = 0 ∧ xα < Kα({}) φ3 ≡ xα = 0 ∧ x β = 1 φ2 1 ≡ xα = 2 ∧ x β = 1 ∧ Kα({α, β}) < 2 φ4 ≡ xα > 0 ∧ x β = 0 φ2 2 ≡ xα = 1 ∧ x β = 1 φ5 ≡ xα = 0 ∧ x β = 0 ∧ Kα({}) = 0. φ2 3 ≡ xα = 1 ∧ x β = 0 ∧ Kα({}) = 0 δ T ⊂ Q T × Q T are
mutually defined by: (q 0 Π , x, ) ∈ Q T and for all q T = (q Π , x, pc) ∈ Q T , for all (q Π , g Π , a P , q Π ) ∈ δ Π such that pc = pc ∧ g Π x = ⊥, then q T = (q Π , a P (x), pc ) ∈ Q T and (q T , q T ) ∈ δ T . If q Π ∈ A Π , then the node is said to be accepting.

For any node (q Π , x, pc) ∈ Q T with q Π = (q P , q B ) ∈ Q P × Q B , pc is the path condition in the form of a constraint over parameters in K. It defines the set of annotated dynamics D that can reach the state (q P , x) from the state in Q P × X associated to the ancestor of (q Π , x, pc). With the process described in the section 3, D itself allows the definition of the set of dynamics of the corresponding GRN which can reach all the states associated to nodes of the node path along the same sequence of traversed states. By construction, path conditions expressed over parameters increase along paths of SET and reduce the number of dynamics compatible with the path under construction.

Biological properties are expressed along infinite sequences, and thus, paths of the product and paths of SET are also infinite. But, by disregarding path conditions, the number of possible nodes in a SET is finite 6 . So, when we are building a new node (q Π , x, pc ) whereas it is descendant of a node (q Π , x, pc) (same vertex in Q Π and same value in X), we stop the analysis of the path; these two nodes are respectively called child node and return node. By construction, the path condition of the child node is included in the path condition of the return node, i.e. all parameter instantiations satisfying the child path condition also satisfy the return path condition.

Thus, by performing a mixed symbolic and numerical execution (parameters in K remain unchanged and state variables in X are evaluated), we can stop the execution procedure of the product Π so that each path of the resulting SET is finite and contains a cycle (starting at the return node and ending with the transition leading to the child node). If there exists an accepting node between the return node and the child node, the path condition is said accepting.

Once all finite SET are built, it remains to compute for which parameter instantiations there exist accepting paths. For that, it suffices to consider every accepting path conditions of the SET associated to Π. Each accepting path condition can be satisfied by (at least) one instantiation of parameters in K, it means that there exists a path in Π going infinitely often through the associated cycle, and thus passing infinitely often by an accepting state. And so there exists a path in the dynamics corresponding to this accepting path condition verifying ¬ϕ.

Thus, instantiations of parameters verifying the conjunction of the negation of every accepting path condition of the SET associated to Π correspond to the dynamics such that there is no path verifying ¬ϕ, in other words, all paths verify ϕ. Note that the obtained dynamics verify ϕ along all paths; if the model must verify ϕ only on at least one path, our approach remains adequate with a small adaptation: to do this, we have to get the disjunction of all accepting path conditions of the SET associated to the product P ⊗ B(ϕ).

Example 8. For the Product in Figure 6, there are two solutions after computation; the corresponding values of parameters are: K α ({}) = 1 or 2, K α ({α}) = 2, K α ({β}) = 0, K α ({α, β}) = 2, K β ({}) = 0 and K β ({α}) = 1. One of the corresponding dynamics (with K α ({}) = 2) is represented in Figure 2.

Algorithm of traversal of SET. Algorithm 1, based on a Depth First Search schema, gives an overview of how we practically compute the accepting path conditions. We use three global variables: the parametric product Π, the list of accepting path conditions acceptingP C, and the list nodesList of SET nodes which have already been analyzed.

Starting with a SET node, line 2 to line 4 test and compute its successors, as explained in the "Symbolic Execution Trees" part of section 4. Three tests are then performed successively. Firstly, if the path condition of the successor node is already known, it cannot provide additional information (the pc becomes more specific every depth call), and we stop the study of this successor (line 5). Secondly, line 7 tests if one of the ancestors of the successor is a return node (ancestor with the same vertex and state). If it is the case, then there is an infinite cycle between them and, if there is an accepting node in that cycle, then the successor node is an accepting return node, and its path condition is added to the list accepting P Cs (lines 8 to 9). Thirdly, if the successor is not a return node, we check (line 10) if the node corresponds to a node in nodesList with the same vertex, the same state and the same or a more general path condition.

else if (q Π , x , pc ) ∈ nodesList with pc ⊂ pc then // no 10 copy node: recall of DFS() DFS((q Π , x , pc ), ancestorsList); 11 nodesList.add((q Π , x , pc ));

12
If it is not the case (no copy node), then the DFS function is recalled with the successor node in argument (line 11), which is then added to nodesList (line 12).

Transient and Stable. Nodes of the tree are of the form (q Π , x, pc) with q Π = (q P , q B ) ∈ Q P × Q B . According to the value of q P (either T or S, from the PGRN), we say that the node is either transient or stable. By construction of the PGRN, the target vertices of all transitions outgoing from a stable vertex are vertices of the same type, hence the appellation stable. Furthermore, the guards of the transitions between the stable vertices are always of the form g Π = g P ∧ g B with g P = , and the assignment of the transitions is a P = id (identity assignment). Thus a specific treatment can be provided for the stable nodes, briefly described in the sequel.

In the algorithm 1, the line 11 can be split in two calls, one for the current function and another for the specific treatment of stable nodes (called if the successor node is stable, called stable root in the sequel). In this case, the second argument of the function, the list of ancestors, is an empty list since none of the previous ancestors (all transient) can be a return node of a stable node. According to the characteristics of the transitions between stable vertices mentioned above, for the treatment of stable nodes there is no need to test if the path condition is already known (line 5, already tested with the corresponding stable root), there is no guard on parameters to symbolically verify and no substitution of levels of expression (lines 3 to 4), and there is no update of the state to do (line 6). Furthermore, if a node is accepting then all the explorations of the SET from the stable root can be stopped; indeed, its path condition is identical to all path conditions of the nodes which can be built from it.

Assessment

The methodology described above has been implemented in a prototype software tool called SPuTNIk. SPuTNIk is written in Java and relies on the Z3 constraint solver [START_REF] De | Z3: An efficient smt solver[END_REF] to check the satisfiability of path conditions during the traversal of SET and on the ltl2ba and LTL2BA4J libraries [START_REF] Gastin | Fast LTL to Büchi automata translation[END_REF][START_REF] Bodden | LTL2BA4J Software[END_REF] to generate a Büchi Automaton of minimal size from an LTL formula. To validate our approach with SPuTNIk, we have considered a common biological case study: the analysis of the genetic network that controls the life cycle of the λ phage virus [START_REF] Thieffry | Dynamical behaviour of biological regulatory networks -II. immunity control in bacteriophage lambda[END_REF]. The λ phage can infect the E. coli bacterium with two different outcomes: either it integrates the genome of the host through a process called lysogeny or it enters a lytic phase where it kills the residing cell to reproduce itself. We based our approach on the λ phage model studied in [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF] by Klarner et al. and composed of four genes, denoted cI, cII, cro and N, and ten interactions described Fig. 7.

cI cro cII N [3, -] [1, -] [ 3 , - ] [2, -] [2, +] [2, -] [2, -] [ 1 , - ] [1, +] [1, +]
Fig. 7: The interaction graph G λ for the λ phage.

Klarner et al. describe biological properties as time series: a sequence of specific states given in the form θ ≡ s 1 , * , s 2 , * , . . . , * , s n where s i is the i th observed state while * denotes a possibly empty sequence of unspecified states. Times series are equivalent to LTL formulas of the form φ ≡ s

1 ∧ F(s 2 ∧ F(• • • ∧ F(s n ) . . . )) (i.
e. only composed of ∧ and F operators). Moreover, the states of time series are fully determined, each level of expression corresponds to a single value. For example for G λ , each state is a quadruple (x cI , x cII , x cro , x N ) ∈ {0, 1, 2} × {0, 1} × {0, 1, 2, 3} × {0, 1}). Given a time series θ and an interaction graph, the goal of Klarner et al. is to find out all models which contain at least one path matching θ (i.e. passing through the states of θ in the correct order). In order to reproduce the same experiment than Klarner et al., we discard the Min/Max constraint for all genes (as it is not supported in [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF]), and we relax the Observation constraint for the specific case where cI is activator of itself (as done in [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF]). We then use SPuTNIk to find out the parameter instantiations corresponding to dynamics which are guaranteed to exhibit either a lytic or a lysogenic phenotype in compliance with series θ 1 and θ 2 (i.e. all models that contain at least one path that satisfies φ 1 and at least one that satisfies φ 2 ).

The obtained results are in accordance: amongst the 7 billions possible models, we obtain the same number (8759) of valid ones as in [START_REF] Klarner | Parameter identification and model ranking of thomas networks[END_REF]. But unlike Klarner et al. our method is not restricted to time series: we can consider any form of LTL formulas and there is no need to fully specify all the levels of expression. For example, it is known that a lytic λ phage can not become lysogenic in the future (and conversely). This knowledge cannot be expressed with time series, but it corresponds to the following LTL formulas: φ 3 ≡ G(lys 2 ⇒ ¬F(lyt 3 )) and φ 4 ≡ G(lyt 3 ⇒ ¬F(lys 2 )). By adding these formulas to the previous, we reduce the number of solutions to 2390.

Conclusion

In this paper we introduced a new methodology for reverse-engineering of genetic network models, based on adaptation of classical LTL model-checking with symbolic execution. In order to find dynamics consistent with biological knowledge, we use the whole extent of LTL to express biological knowledge in terms of constraints over time. Instead of checking each dynamics of the GRN, we propose a method which performs checking with a novel formalism, the Parametric GRN, a compact (symbolic) representation of all the dynamics associated to an interaction graph within a single structure. From the Parametric GRN and LTL formulas, our algorithm processes parameters, defining the dynamics, as symbols in order to avoid combinatorial explosion. The solutions are in the form of a set of constraints that the parameters must fulfill. Such analysis has been carried out through the SPuTNIk tool, a prototype software of the proposed method. We are working on a parallel version of SPuTNIk, based on the splitting of the Parametric Product into strongly connected components in order to detect the accepting cycles in each component.

  Klarner et al. distinguish the following states: init ≡ [0000], lyt 1 ≡ [0021], lyt 2 ≡ [0020], lyt 3 ≡ [0030], lys 1 ≡ [2101] and lys 2 ≡ [2000], belonging to time series θ 1 ≡ init, * , lyt 1 , * , lyt 2 , * , lyt 3 and θ 2 ≡ init, * , lys 1 , * , lys 2 , which correspond to evolution towards lytic and lysogenic phases. The equivalent LTL counterparts for θ 1 and θ 2 are respectively φ 1 ≡ init ∧ F(lyt 1 ∧ F(lyt 2 ∧ F(lyt 3 ∧ F(lyt 2 )))) and φ 2 ≡ init ∧ F(lys 1 ∧ F(lys 2 )).

A transition system (E, R) verifies that R is a binary relation on E × E.

The number of possible parameters instantiation is equal to g∈G (mg + 1) 2 |G -(g)| .

it is bounded by the product of all combinations of levels of expression, the number of vertices of the Büchi automaton and the number of vertices (2) of the PGRN.

Algorithm 1: Overview of DFS((q Π , x, pc), ancestorsList)


Astrid Lièvre

Bérèngere Ouine

Jim Canet

Aurélie Cartier

Yael Amar

Wulfran Cacheux

Odette Mariani

Rosine Guimbaud

Janick Selves

Thierry Lecomte

Berengere Ouine

Serge Guyetant

Ivan Bieche

Frédérique Berger

Leanne De Koning
email: leanne.de-koning@curie.fr

Protein biomarkers predictive for response to anti-EGFR treatment in RAS wild-type metastatic colorectal carcinoma

Keywords: colorectal cancer, resistance, EGFR, biomarkers, cell signaling, RPPA, PI3K

Introduction

Colorectal cancer (CRC) is the second leading cause of cancer-related deaths, with 1.4 million new cases world-wide in 2012. The prognosis of CRC is mainly related to the presence of metastasis: around 25% of patients present with metastasis upon diagnosis and around 50% of patients that are treated for localized CRC will develop metastases during the course of disease. Despite the advances in early diagnosis and treatment achieved in the past 20 years, prognosis of metastatic CRC (mCRC) remains relatively poor, with a 5-year relative survival rate of about 12% (American Cancer Society).

Epidermal Growth Factor Receptor (EGFR or ERBB1) is a cell membrane receptor that belongs to the family of receptor tyrosine kinases [START_REF] Arteaga | ERBB receptors: from oncogene discovery to basic science to mechanism-based cancer therapeutics[END_REF]. Upon binding of various ligands such as EGF, the receptor is activated and induces the activation of downstream signaling pathways, including PI3K/AKT, MEK/ERK, Jak/Stat and JUNK pathways, which contribute to tumorigenesis. Overexpression or activating mutations of EGFR occur in many cancer types, among which CRC. The development of monoclonal antibodies directed against EGFR, such as cetuximab and panitumumab, has significantly improved CRC outcome, both in the context of chemoresistant tumors [START_REF] Amado | Wild-type KRAS is required for panitumumab efficacy in patients with metastatic colorectal cancer[END_REF][START_REF] Karapetis | K-ras mutations and benefit from cetuximab in advanced colorectal cancer[END_REF]) and as a first-line treatment [START_REF] Bokemeyer | Fluorouracil, leucovorin, and oxaliplatin with and without cetuximab in the first-line treatment of metastatic colorectal cancer[END_REF][START_REF] Van Cutsem | Cetuximab and chemotherapy as initial treatment for metastatic colorectal cancer[END_REF]). However, only patients having a tumor without mutations in KRAS and NRAS can benefit from anti-EGFR therapy [START_REF] Douillard | RAS mutations in colorectal cancer[END_REF]. Indeed, KRAS and NRAS operate downstream of EGFR in the RAS/MAPK signaling pathway and their mutation activates the pathway independently of EGFR status. KRAS and NRAS mutations are frequent, occurring in around 50% of CRCs, and their sequencing is therefore mandatory before administration of anti-EGFR treatment.

However, even among the patients with wildtype KRAS and NRAS, only 20 to 30 % respond to the anti-EGFR treatment monotherapy [START_REF] Price | Final results and outcomes by prior bevacizumab exposure, skin toxicity, and hypomagnesaemia from ASPECCT: randomized phase 3 non-inferiority study of panitumumab versus cetuximab in chemorefractory wild-type KRAS exon 2 metastatic colorectal cancer[END_REF]) and 65-70% to anti-EGFR combined with chemotherapy [START_REF] Heinemann | A studylevel meta-analysis of efficacy data from head-to-head first-line trials of epidermal growth factor receptor inhibitors versus bevacizumab in patients with RAS wild-type metastatic colorectal cancer[END_REF], suggesting that other molecular mechanisms of resistance exist. The identification of additional markers of resistance would allow to better select those patients that could benefit from anti-EGFR therapy and avoid inefficient and potentially toxic treatment of the other patients. CRC cell lines, xenografts and, less frequently, patient samples have been used to tackle this question. Multiple studies have shown that activation of the signaling pathways downstream of EGFR, induced by genetic alterations such as PTEN deletions, PIK3CA mutations or MET activation, constitute an important mechanism of primary and acquired resistance towards anti-EGFR [START_REF] Bardelli | Molecular mechanisms of resistance to cetuximab and panitumumab in colorectal cancer[END_REF][START_REF] Troiani | Increased TGF-alpha as a mechanism of acquired resistance to the anti-EGFR inhibitor cetuximab through EGFR-MET interaction and activation of MET signaling in colon cancer cells[END_REF][START_REF] Bajpe | A chromatin modifier genetic screen identifies SIRT2 as a modulator of response to targeted therapies through the regulation of MEK kinase activity[END_REF][START_REF] Luraghi | MET signaling in colon cancer stem-like cells blunts the therapeutic response to EGFR inhibitors[END_REF][START_REF] Song | Cetuximab-induced MET activation acts as a novel resistance mechanism in colon cancer cells[END_REF][START_REF] Van Emburgh | Acquired resistance to EGFR-targeted therapies in colorectal cancer[END_REF].

HER2 amplification or mutation has also been associated with anti-EGFR resistance in CRC xenografts [START_REF] Bertotti | A molecularly annotated platform of patient-derived xenografts ("xenopatients") identifies HER2 as an effective therapeutic target in cetuximab-resistant colorectal cancer[END_REF][START_REF] Yonesaka | Activation of ERBB2 signaling causes resistance to the EGFR-directed therapeutic antibody cetuximab[END_REF][START_REF] Bertotti | The genomic landscape of response to EGFR blockade in colorectal cancer[END_REF]. Furthermore, amplifications or mutations of FGFR1, PDGFRA, MAP2K1 have been described [START_REF] Bertotti | The genomic landscape of response to EGFR blockade in colorectal cancer[END_REF], as well as the deregulation of several microRNAs. Yet, besides RAS, no other marker of resistance has been validated so far for clinical practice. In addition, recent data suggest that right-sided mCRC is more resistant to anti-EGFR therapy than left-sided mCRC [START_REF] Moretto | Location of Primary Tumor and Benefit From Anti-Epidermal Growth Factor Receptor Monoclonal Antibodies in Patients With RAS and BRAF Wild-Type Metastatic Colorectal Cancer[END_REF], Tejpar et al. 2016[START_REF] Holch | The relevance of primary tumour location in patients with metastatic colorectal cancer: A meta-analysis of first-line clinical trials[END_REF], but the biology underlying this difference remains elusive.

We hypothesize that many different genetic or epigenetic alterations can be involved in anti-EGFR resistance and that the crucial end point resides in the activation of downstream signaling pathways. The activation of these pathways would thus be a better and more universal predictor of resistance than each genetic alteration separately. However, large-scale protein data on CRC patient samples with clinical follow-up is currently missing. For this reason, we here decided to analyze the expression and the activation of a large panel of proteins involved in cell signaling pathways, using Reverse Phase Protein Arrays (RPPA) on a multi-center French cohort of RAS wild-type mCRCs, both left-and right-sided, treated with anti-EGFR treatment.

Methods and materials

Patient samples

Patients (n=53) with metastatic chemoresistant CRC were treated with anti-EGFR therapy (cetuximab or panitumumab), alone or in association with chemotherapy, at Institut Curie (Paris, France), CHU of Toulouse (France) or CHRU of Tours (France). Patients could be included in the study if tumor response to anti-EGFR could be specifically evaluated i.e. patients treated with a combination of anti-EGFR and chemotherapy who previously progressed on the same chemotherapy component (including those who had progressed on an oxaliplatin-based adjuvant chemotherapy then on a first-line irinotecan-based chemotherapy and those who had progressed on a first-line FOLFIRINOX tritherapy), or patients treated with a monotherapy of anti-EGFR or patients who initially progressed on a first-line combination of anti-EGFR and chemotherapy. According to French regulations, patients were informed of research performed with the biological specimens obtained during their treatment and did not express opposition. This retrospective study was reviewed and approved by the Ethics Committee of the Institut Curie. Time from sample resection to sample freezing was less than 30 minutes in most cases, and always less than 60 minutes. Samples were stored in secured -80°C freezers. For this retrospective study, four 50µm thick frozen tissue sections of the primary tumor, obtained before administration of anti-EGFR treatment and containing at least 50% of tumor cells, were sent to the RPPA platform of Curie.

Sequencing of KRAS, NRAS, BRAF and PIK3CA was performed independently on 3 of the 28 regional molecular genetics platforms constituting the national network of public laboratories dedicated to molecular oncology tests in France that has been certified by the French National Cancer Institute (INCa).

According to the INCa quality assurance program, these platforms have used one of the recommended sequencing techniques with detection sensitivity between 5% and 10% of mutated cells (allelic hybridization using HRM followed by Sanger Sequencing or by pyrosequencing). From the 53 samples, 7 were removed due to low RPPA signals, probably reflecting protein degradation. 12 more samples had to be excluded from further analysis, because they did not comply with inclusions criteria (6 had a KRAS mutation that was not initially determined since these patients were diagnosed before KRAS sequencing became a prerequisite for anti-EGFR treatment, for 2 tumors cellularity was too low, 2 patients were responders to a first-line combination of anti-EGFR and chemotherapy so that we could not determine the specific response to the anti-EGFR, and 2 tumors were not CRC). Thus, 34 samples were kept for further analysis.

Reverse Phase Protein Arrays

Samples were disrupted in Laemmli buffer (50 mM Tris pH =6.8, 2% SDS, 5% glycerol, 2 mM DTT, 2.5 mM EDTA, 2.5 mM EGTA, 1x HALT Phosphatase inhibitor (Perbio 78420), Protease inhibitor cocktail complete MINI EDTA-free (Roche 1836170, 1 tablet/10 mL), 2 mM Na3VO4 and 10 mM NaF), using a Precellys (Bertin). Extracts were then boiled for 10 min at 100°C, sonicated to reduce viscosity and centrifuged 10 min at 15000 rpm. The supernatant was harvested and stored at -80°C. Protein concentration was determined (Pierce BCA reducing agent compatible kit, ref 23252). Samples were printed onto nitrocellulose covered slides (Supernova, Grace Biolabs) using a dedicated arrayer (2470 arrayer, Aushon Biosystems). Five serial dilutions, ranging from 1500 to 94 µg/ml, and three technical replicates per dilution were printed for each sample. Arrays were labeled with 86 specific antibodies (see supplementary Table 1 for a complete list of antibody references) or without primary antibody (negative control), using an Autostainer Plus (Dako). Briefly, slides were incubated with avidin, biotin and peroxidase blocking reagents (Dako) before saturation with TBS containing 0.1% Tween-20 and 5% BSA (TBST-BSA). Slides were then probed overnight at 4°C with primary antibodies diluted in TBST-BSA. After washes with TBST, arrays were probed with horseradish peroxidase-coupled secondary antibodies (Jackson ImmunoResearch Laboratories) diluted in TBST-BSA for 1h at RT. To amplify the signal, slides were incubated with Bio-Rad Amplification Reagent for 15 min at RT. The arrays were washed with TBST, probed with Alexa647-Streptavidin (Molecular Probes) diluted in TBST-BSA for 1h at RT and washed again in TBST. For staining of total protein, arrays were incubated 15 min in 7% acetic acid and 10% methanol, rinsed twice in water, incubated 10 min in Sypro Ruby (Invitrogen) and rinsed again. The processed slides were dried by centrifugation and scanned using a GenePix 4000B microarray scanner (Molecular Devices). Spot intensity was determined with MicroVigene software (VigeneTech Inc). All primary antibodies used in RPPA have been previously tested by Western Blotting to assess their specificity for the protein of interest.

Samples with low overall signal or with aberrant dilution curves, which are often indicative of protein degradation, were discarded during quality control. For each antibody, the median signal intensity was at least three times higher than the background array without primary antibody. Raw data were normalized using Normacurve [START_REF] Troncale | NormaCurve: a SuperCurve-based method that simultaneously quantifies and normalizes reverse phase protein array data[END_REF], which normalizes spot-wise for a negative control slide (labeled without any primary antibody, provides the fluorescent background) and a slide labeled with a total protein stain (serves as a loading control). Next, Normacurve draws the antibody response curve for each array. Each sample, including the five serial dilutions and the replicates, are aligned onto this curve to calculate one normalized value per sample. These normalized values are then used for statistical analysis.

Statistical Analysis

Response to anti-EGFR treatment was determined using RECIST criteria [START_REF] Therasse | New guidelines to evaluate the response to treatment in solid tumors. European Organization for Research and Treatment of Cancer, National Cancer Institute of the United States, National Cancer Institute of Canada[END_REF].

Continuous variables are described as mean and standard deviations and qualitative data are presented as a number and percentage of sample size. The association between clinical variables and proteins was determined by Wilcoxon or Kruskal-Wallis tests. Association between discrete variables was tested using Chi-square or Fisher's exact test. Hierarchical clustering was performed using Ward metrics and Pearson correlation and represented as a heatmap.

Protein expression was divided in high expression and low expression with the cut-off at the median expression level.

The RECIST Criteria were restricted to a dichotomous output, where complete or partial response were considered as a response and stabilization or progression as a non-response. Univariate logistic regression analysis was performed and odds ratios were calculated to measure the association between the expression level of a protein and the chance to respond to treatment. A multivariate logistic regression analysis was performed including significant proteins (p<0.05) with a stepwise procedure.

Overall survival is defined as the time between diagnosis of the metastasis and the date of death.

Patients still alive at the moment of analysis were censured at the date of last follow-up. Survival curves were estimated using Kaplan-Meier and compared with Log-rank tests.

Univariate and multivariate Cox proportional hazard models were performed to determine the variables that impacts OS. Only variables with a significant p value (p<0.05) were included in a multivariate stepwise procedure using the Cox model. P-values below 0.05 were considered significant. All the analyses were performed using R software version 3.3.0 (R Core Team 2015).

Results

Patient characteristics and response to anti-EGFR therapy according to primary tumor location

Final RPPA data was obtained for 34 CRCs and 86 antibodies. The antibodies were selected according to signaling pathways that have been put forward as being involved in EGFR signaling and anti-EGFR resistance (various RTKs, PI3K/Akt pathway, MAPK pathways) or more generally in CRC carcinogenesis (Wnt/Notch) and chemotherapy response (apoptosis, cell proliferation, DNA repair, angiogenesis).

Clinical characteristics are summarized in table 1 and antibodies are listed in Supplementary Table 1.

Median follow-up of these 34 patients was 15.5 months from the beginning of anti-EGFR treatment (range 3 -45.6 months). All patients received anti-EGFR as second or third line treatment, except for four patients who received anti-EGFR in combination with chemotherapy as a first line treatment but did not respond. All patients of this cohort were thus shown to be chemoresistant, and any observed response to therapy could be attributed solely to the addition of anti-EGFR treatment. 11 patients showed a partial response, 11 patients showed stabilization and 12 patients showed disease progression.

No complete response was observed. Among the 11 patients with stabilization, the median time to progression was 8 months (range: 3-13 month) versus 16 months (range: 5-46) in patients with partial response and 2 months (range:1-3 months) in patients with progressive disease respectively. Given the small sample size of patients with stabilization, we did not perform a further categorization of stabilized patients according to duration of stabilization, which would limit the significance of the results. Five patients showed a PIK3CA hot spot activating mutation (3 patients have p.E545K, one patient p.E542K and one patient p.N1044K) and one other patient a BRAF V600E mutation. Biopsies were obtained before initiation of the anti-EGFR therapy.

We analyzed if the response to treatment, measured according to RECIST criteria, was associated with any of the clinical parameters. The RECIST score was not significantly different according to sex (men versus female), age (<50 versus ≥50) or the molecule of anti-EGFR treatment (cetuximab versus panitumumab) that was administered. As expected, the absence of objective response was associated to worse overall survival (data not shown). Interestingly, and in accordance with recent findings [START_REF] Moretto | Location of Primary Tumor and Benefit From Anti-Epidermal Growth Factor Receptor Monoclonal Antibodies in Patients With RAS and BRAF Wild-Type Metastatic Colorectal Cancer[END_REF], Tejpar et al. 2016[START_REF] Holch | The relevance of primary tumour location in patients with metastatic colorectal cancer: A meta-analysis of first-line clinical trials[END_REF]) the response to anti-EGFR treatment was better in leftsided CRC (descending colon and rectum), where 46% of patients show partial response, than in rightsided CRC (ascending colon), where none of the patients shows partial response (p= 0.03).

Protein expression according to patient and tumor characteristics

We studied if certain proteins, measured by RPPA, were associated with clinical parameters. None of the measured proteins showed a significant association with sex, age or the number of metastatic sites (1 vs ≥2). Next, we compared the protein profiles of left-sided versus right-sided CRC. Out of 86 protein analyzed, 76 proteins are differentially expressed or activated between the two locations (p<0.05), demonstrating an important difference in terms of signaling pathway activation. The phosphorylated proteins that are differentially expressed between left and right sided CRC are enriched in HER2/HER3 signaling (p=0.014) and in the Wnt pathway (p=0.018) (figure 1A), which both seem to be more activated in left sided CRC as previously suggested [START_REF] Kim | Sex-and gender-specific disparities in colorectal cancer risk[END_REF] (figure 1B and1C).

Protein biomarkers predictive for response to anti-EGFR therapy

Using the RECIST criteria, we searched for (phospho-)proteins associated with response to anti-EGFR antibodies that could thus constitute potential predictive biomarkers. Because of small group sizes and to better identify biomarkers the most predictive of response to anti-EGFR, we compared partial response (n=11) versus stabilization + progression (n=23). Higher levels of Phospho-Akt (Ser473) (p=0.01), HER2 (p=0.03), PKC delta (p=0.03), Phospho-HER4 (p=0.04), Phospho-p70S6kinase (p=0.05) and 4EBP1 (p=0.05) are associated with response to treatment in univariate analyses (Figure 2). In addition, several proteins show a trend towards significance: higher levels of phospho-EGFR (p=0.06), FGFR4 (p=0.06), phospho-GSK3 (p=0.07) and p53 (p=0.07) are associated with a better response to treatment (Figure 2). Multivariate analyses allow evidencing Phospho-Akt (Ser473) as the dominating biomarker for response (OR=5.5, CI95% [1.6 ;34.2]).

Prognostic factors associated with survival

Next, we studied which characteristics were associated with overall survival, defined as the time between the diagnosis of the metastasis and death. Using unsupervised hierarchical clustering, the 86 analyzed proteins allow a separation of the 34 CRCs into two clearly distinguishable groups (Supplementary Figure 1). The two clusters do not separate the patients according to response, location of the primary tumor (left-or right-sided, or rectum versus colon) or the center of origin of the samples.

In addition, the two clusters do not display a significant difference in survival (data not shown). Leftsided CRC seem to have a better overall survival compared to right-sided CRC (survival rate at 24 months of 78% (CI95% [63% ; 98%]) and 36% (CI95%[13% ; 99%]) respectively), although the log rank test is not significant (Supplementary Figure 2). Next, we studied which individual (phospho-)proteins are associated with overall survival. We found that a high expression of EGFR (p=0.01), Phospho-EGFR (Tyr1173) (p=0.03) and HER3 (p=0.03) are associated with a better survival (Figure 3). The expression of Phospho-EGFR and HER3 proteins is highly correlated (correlation coefficient R= +0.82, p<0.001) and they show identical survival curves (Figure 3B and 3C).

Both proteins also correlate with total EGFR (R= 0.62 and R= 0.61 respectively; p<0.001). Higher expression of FGFR3 (p=0.05), phospho-4EBP1 (p=0.05), p53 (p=0.06) and phospho-HER3 (p=0.06) show a trend towards an association with better survival, without reaching significance.

Discussion

Only a fraction of mCRC patients having wild-type KRAS and NRAS benefit from anti-EGFR treatment, suggesting the presence of additional molecular characteristics leading to primary resistance. Many genetic alterations have been shown to induce resistance in cell lines or xenograft models. Here, we hypothesized that all these alterations will ultimately lead to the activation of cell signaling pathways that can be measured at the protein level.

We performed the largest targeted proteomics study published so far in terms of analyzed proteins on a small but well annotated cohort of 34 KRAS and NRAS WT mCRC samples from patients that received anti-EGFR therapy, with the aim to identify predictive markers of sensitivity or resistance.

We identified several (phospho-)proteins that are predictive for response to treatment or for overall survival in mCRC patients. Although the observed differences per protein are small and not highly significant due to the small study size, the identified proteins reveal interesting patterns. Indeed, we show that, independently of the line of treatment, patients with higher expression and activation of EGFR and HER3 membrane receptors have a better overall survival. HER3, which lacks a functional kinase domain, heterodimerizes with EGFR or with HER2 to produce a potent signaling complex [START_REF] Jura | Structural analysis of the catalytically inactive kinase domain of the human EGF receptor 3[END_REF]. Targeting EGFR and HER3 concomitantly is a current lead in CRC [START_REF] Juric | Safety and Pharmacokinetics/Pharmacodynamics of the First-in-Class Dual Action HER3/EGFR Antibody MEHD7945A in Locally Advanced or Metastatic Epithelial Tumors[END_REF][START_REF] Temraz | Dual targeting of HER3 and EGFR in colorectal tumors might overcome anti-EGFR resistance[END_REF]. Thus, actively signaling EGFR is associated with a better overall survival in these patients receiving anti-EGFR treatment, probably because these tumors are more dependent on EGFR signaling and thus more sensitive to its inhibition.

If we look more specifically at the response to treatment, as defined by the RECIST criteria, a broader picture appears. The proteins that associate with a better response to therapy are mostly associated with an activation of tyrosine kinases (EGFR, HER2, HER4 and FGFR4) and the downstream signaling through the PI3K/Akt/mTOR pathway (GSK3 and 4EBP1 which are downstream of Akt, p70S6K which is downstream of mTOR, PKCɗ which can be activated by EGFR), thus revealing a complex regulatory network (Figure 4). Thus, patients with activated cell surface receptors (notably EGFR) and PI3K pathway are more likely to respond to anti-EGFR therapy. Our data confirm and extent previously reported observations on the predictive value of Phospho-EGFR and Phospho-Akt [START_REF] Van Schaeybroeck | Epidermal growth factor receptor activity determines response of colorectal cancer cells to gefitinib alone and in combination with chemotherapy[END_REF][START_REF] Harle | Expression of pEGFR and pAKT as response-predictive biomarkers for RAS wild-type patients to anti-EGFR monoclonal antibodies in metastatic colorectal cancers[END_REF].

In CRC carrying a mutation in KRAS or NRAS, this mutation confers resistance to anti-EGFR treatment by activating MAPK and PI3K pathways. In our patient cohort of RAS WT CRCs, we observe active EGFR signaling and downstream PI3K pathway activation. Interestingly, we do not identify components of the MAPK pathway as being predictive for response to treatment, suggesting that in RAS WT patients the PI3K pathway is the predominant pathway that explains variability in response to anti-EGFR therapies.

Several potential mechanisms could be at the origin of the EGFR activation in these RAS WT tumors. First, the overexpression of EGFR ligands and notably Epiregulin and Amphiregulin, which activate EGFR, has been associated with a better response to anti-EGFR therapy in RAS WT CRC [START_REF] Khambata-Ford | Expression of epiregulin and amphiregulin and K-ras mutation status predict disease control in metastatic colorectal cancer patients treated with cetuximab[END_REF][START_REF] Baker | Tumour gene expression predicts response to cetuximab in patients with KRAS wild-type metastatic colorectal cancer[END_REF][START_REF] Jonker | Epiregulin gene expression as a biomarker of benefit from cetuximab in the treatment of advanced colorectal cancer[END_REF][START_REF] Seligmann | Combined Epiregulin and Amphiregulin Expression Levels as a Predictive Biomarker for Panitumumab Therapy Benefit or Lack of Benefit in Patients With RAS Wild-Type Advanced Colorectal Cancer[END_REF]). Second, mutations in genes such as PIK3CA, PTEN, EGFR, and ERBB2, were recently found predictive for anti-EGFR therapy in 31% of RAS WT tumors [START_REF] Rankin | Broad Detection of Alterations Predicted to Confer Lack of Benefit From EGFR Antibodies or Sensitivity to Targeted Therapy in Advanced Colorectal Cancer[END_REF]). In our cohort, mutation status was determined for BRAF and PIK3CA. A single patient had a BRAF-mutated tumor and had progressive disease as expected. Indeed, the V600E mutation in the gene that encodes BRAF, which acts downstream of RAS, is known to confer resistance to anti-EGFR therapy and a very poor prognosis of CRC patients [START_REF] Pietrantonio | Predictive role of BRAF mutations in patients with advanced colorectal cancer receiving cetuximab and panitumumab: a meta-analysis[END_REF][START_REF] Rowland | Meta-analysis of BRAF mutation as a predictive biomarker of benefit from anti-EGFR monoclonal antibody therapy for RAS wild-type metastatic colorectal cancer[END_REF]. PIK3CA mutations were previously also suggested to lead to resistance to anti-EGFR therapy [START_REF] Sartore-Bianchi | PIK3CA mutations in colorectal cancer are associated with clinical resistance to EGFR-targeted monoclonal antibodies[END_REF], although large scale meta-analyses suggest that this is true only for mutations in exon 20 and not for mutations in exon 9 [START_REF] De Roock | Effects of KRAS, BRAF, NRAS, and PIK3CA mutations on the efficacy of cetuximab plus chemotherapy in chemotherapy-refractory metastatic colorectal cancer: a retrospective consortium analysis[END_REF]. We detected activating PIK3CA hot spot mutations in 5 patients: 4 mutations in exon 9 (three p.E545K and one p.E542K mutation) and one in exon 20 (p.N1044K). From these 5 patients, three had a partial response (including the patient with the exon 20 mutated tumor), one had a stable disease and one a progressive disease.

Although the numbers are too small for statistical analysis, there is thus no indication in our cohort that PIK3CA mutation associates with anti-EGFR resistance. On the contrary, PIK3CA mutations could be one method to activate the PI3K pathway, which we find associated with a better response to therapy.

During the course of our study, large scale meta analyses revealed that left-sided and right-sided CRC do not respond similarly to anti-EGFR therapy [START_REF] Moretto | Location of Primary Tumor and Benefit From Anti-Epidermal Growth Factor Receptor Monoclonal Antibodies in Patients With RAS and BRAF Wild-Type Metastatic Colorectal Cancer[END_REF], Tejpar et al. 2016[START_REF] Holch | The relevance of primary tumour location in patients with metastatic colorectal cancer: A meta-analysis of first-line clinical trials[END_REF].

We here confirm that left-sided CRC has a better response to anti-EGFR therapy, a tendency towards better survival and a very different profile of protein expression with notably more Wnt and ErbB signaling activation.

In conclusion, we identified activated EGFR and HER3 as biomarkers predictive for a better overall survival in patients treated by anti-EGFR therapy. Response to treatment, on the other hand, was associated with several markers that converge to active EGFR signaling and in particular the PI3K pathway. Validation of these markers by immunohistochemistry on a large panel of samples would therefore be a crucial step forwards to improved patient stratification and personalized medicine in RAS wildtype CRC. Titles and Legends to Figure

Figure 1 :

 1 Figure 1: Differences in pathway activation between left-sided versus right-sided CRC. A. Phosphorylated

Figure 2 :

 2 Figure 2: Distribution of (phospho)proteins that are differentially expressed according to the response to

Figure 3 :

 3 Figure 3: Kaplan Meier curves of overall survival according to expression levels of EGFR (A), Phospho-

Figure 4 :

 4 Figure 4: Simplified scheme of signaling interaction network between the proteins that are associated

Table 1 : Description of Clinical Variables of mCRC samples

 1

	N	%

Acknowledgements

We thank Stéphane Liva and Patrick Poullet for bioinformatics support. This work was financed by the Institut Curie Evaluation Committee of Translational Research (CEST) and Institut Carnot Curie-Cancer 2013.

Correlation between mutational data, response to treatment and protein expression

In our cohort, a single patient had a BRAF V600E mutation and had progressive disease. Five patients had PIK3CA mutations. From these 5 patients, three had a partial response; one had a stable disease and one a progressive disease. Although the numbers are too small for statistical analysis, there is thus no indication in our cohort that PIK3CA mutation associates with anti-EGFR resistance. Despite the lack of statistical power, we addressed if the 5 PIK3CA mutated tumors indeed show a higher activation of the Akt pathway than PIK3CA WT tumors. The ratios of Phospho-Akt (Thr308)/Akt, of phospho-p70S6K/p70S6K and of P-PKCα/PKCα (a downstream target of PDK1 and mTOR) were slightly increased when PIK3CA is mutated (p=0.11; p=0.14 and p=0.10, respectively), without reaching significance though (Supplementary Figure 3).

Conflict of Interest statement

The authors have no conflicts of interest to declare.

Supplementary data

Supplementary Figure 1: Hierarchical clustering of the RPPA data shows two distinct groups of samples that are not related to response to therapy (color bar 1), to whether the tumor is right-or left-sided (color bar 2), to the center of origin (color bar 3) or to the type of tumor (color bar 4). Samples are horizontal, antibodies are vertical. White: missing data.

Supplementary