Fabian Gruber

Manuel Selva

Diogo Sampaio

Christophe Guillon

Louis-Noël Pouchet

Fabrice Rastello

Building of a Polyhedral Representation from an Instrumented Execution: Making Dynamic Analyses of non-Ane Programs Scalable

Keywords: Dynamic analysis, Polyhedral representation, Loop optimization, Polyhedral optimization, Trace compression, Création Analyse dynamique, Représentation polyédrique, Optimisation de boucle, Optimisation polyédrique, Compression de trace

The polyhedral model has been successfully used in production compilers. Nevertheless, only a very restricted class of applications can benet from it. Recent proposals investigated how runtime information could be used to apply polyhedral optimization on applications that do not statically t the model. In this work, we go one step further in that direction. We propose a dynamic analysis that builds a compact polyhedral representation from a program execution. It is able to accurately detect ane dependencies and xed-stride memory accesses in programs. The analysis scales to real-life applications, which often include some non-ane dependencies and accesses in otherwise ane code. This is enabled by a safe ne-grain polyhedral over-approximation mechanism applied to each analyzed expression. We evaluate our analysis on the entire Rodinia benchmark suite, enabling accurate feedback about potential for complex polyhedral transformations.

Introduction

The most eective program transformations for improving performance or energy consumption are typically based on rescheduling of instructions so as to expose data locality and/or parallelism. The two main challenges for those kinds of optimizations are what transformations may be applied without breaking the program semantics and where in the program should the optimizations be applied so as to maximize the impact on the overall program performance. The polyhedral model [START_REF] Feautrier | Polyhedron model[END_REF], usually applied to static compilation, is a powerful tool for nding and applying such rescheduling optimizations. Compilers using the polyhedral model [START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF][START_REF] Grosser | Polly -performing polyhedral pptimizations on a low-level intermediate representation[END_REF] leverage precise information about data and control-ow dependencies to determine a sequence of loop transformations. These loop transformations aim to improve temporal-and spatial locality and uncover both coarse (i.e., thread) and ne-grain (i.e., SIMD) parallelism.

In practice, the static analysis used to recover the dependence information required for these transformations can only be applied to programs written in a very restrictive style with no function calls, only arrays as data structures, only very simple conditional statements and no indirections [START_REF] Collard | Fuzzy array dataow analysis[END_REF]. Dynamic analysis frameworks [START_REF] Faxén | Embla -data dependence proling for parallel programming[END_REF] address this limitation by reasoning on a particular execution of the program. The feedback provided by existing frameworks mainly informs about the absence of dependencies along some loop in the original program, highlighting opportunities for parallelism [START_REF] Ketterlin | Proling data-dependence to assist parallelization: Framework, scope, and optimization[END_REF][START_REF] Wang | Integrating proledriven parallelism detection and machine-learning-based mapping[END_REF][START_REF] Vandierendonck | The paralax infrastructure: automatic parallelization with a helping hand[END_REF][START_REF] Tournavitis | Semi-automatic extraction and exploitation of hierarchical pipeline parallelism using proling information[END_REF] or SIMD vectorization [START_REF] Holewinski | Dynamic trace-based analysis of vectorization potential of applications[END_REF].

In this work, we propose a dynamic analysis that does not merely prove the absence of dependencies, but recovers their structure in a form that ts the polyhedral model. We call this analysis the folding-based analysis. It can accurately detect polyhedral dependencies in programs.

The analysis scales to real-life applications, which often include some non-ane dependencies in otherwise ane code. For that, we propose a safe ne-grain polyhedral over-approximation mechanism for these dependencies. That is, our analysis emits a compact program representation allowing a classic polyhedral compiler to nd a wide range of possible transformations. Our analysis also allows detecting the presence of xed-stride memory accesses. This information is useful for exposing potential for vectorization and loop transformations to improve spatial locality. The contribution of this paper is the folding-based analysis which: builds a compact polyhedral program representation from a program execution, enabling polyhedral compilers to be applied, that is, to provide feedback about the potential of complex polyhedral transformations [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF]; captures information useful for polyhedral optimizers such as properties of dependencies, data ow, memory accesses, and scalar evolution in a uniform manner; and scales to real-life applications by widening of non-ane expressions with a safe polyhedral over-approximation. This paper is organized as follows. Section 2 illustrates the context of our work through a case study. It then continues with an in-depth description of the interface in Section 3, followed by the core algorithm used by our analysis in Section 4. Section 5 evaluates our approach by applying it to the entire Rodinia [START_REF] Che | Rodinia: A benchmark suite for heterogeneous computing[END_REF][START_REF] Che | A characterization of the rodinia benchmark suite with comparison to contemporary cmp workloads[END_REF] benchmark suite. Section 6 discusses related work.

Finally, Section 7 concludes the paper and provides future perspectives.

Motivational Scenario

This section introduces the problem tackled by this work using a concrete example. For this, we use backprop, a benchmark from the Rodinia benchmark suite [START_REF] Che | Rodinia: A benchmark suite for heterogeneous computing[END_REF][START_REF] Che | A characterization of the rodinia benchmark suite with comparison to contemporary cmp workloads[END_REF]. backprop is a supervised RR n°9244 for (j = 1; j <= n2; j++) { // For each unit in second layer float sum = 0.0; // Compute weighted sum of its inputs for (k = 0; k <= n1; k++) sum += conn[k][j] * l1[k]; l2[j] = squash(sum); } Figure 1: A compute intensive kernel in backprop learning method used to train articial neural networks. We choose this benchmark because it is relatively short and simple, and still has behavior that cannot be handled accurately by existing analysis. For illustrative purposes, we focus on the compute kernel shown in Fig. 1. This kernel is also used as a running example throughout the rest of the paper.

As detailed below, while the source code appears ane, in fact taken in context of the full program, a compiler would assume overly pessimistic may-alias dependences due to how arrays are allocated. Dynamic analysis is required here to perform accurate polyhedral optimization on this kernel. To that end, we have developed a complete proling tool-chain [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF]. This tool-chain works on compiled binaries and provides suggestions for loop transformations. It is made of three parts. The front-end, which instruments the binary to get information from an execution of the program. The folding-based analysis, which consumes this information in a streaming fashion to build a compact polyhedral program representation along with information on the behavior of memory accesses. And nally, the back-end, which uses this representation to nd interesting loop transformations. This paper focuses on the folding-based analysis, while the front-end and back-end have been presented in prior work [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF].

2.1

Example problem: backprop Back to Fig. 1, at rst glance, all loop bounds and memory addresses seem to be ane expressions of loop invariant parameters, n2 and n1, and enclosing loop iterators, j and k.

This kernel should thus be a perfect target for polyhedral compilers such as LLVM-Polly [START_REF] Grosser | Polly -performing polyhedral pptimizations on a low-level intermediate representation[END_REF] or GCC-Graphite [START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF]. Nevertheless, this is not the case because the conn object is not a two-dimensional array, but an array of pointers, each allocated by a separate call to malloc as illustrated in Fig. 2. Because of that, the compiler has to generate two load instructions for the access conn[k][j] as shown in Fig. 4. The rst one loads conn[k] into a temporary variable tmp and the second one loads the value of tmp[j] into another temporary used in the multiplication. In other words, while one would expect a static compiler to see a single ane access to a two-dimensional array, there are actually two indirect accesses. Another problem is that it is impossible to know statically if the pointers conn[k], l1, or l2 alias, that is, whether they refer to the same object at runtime. Due to this, a static compiler has to conservatively assume there is a dependence between the write on Line 5 and the reads on Line 4. This dependence prevents any transformation of the outer loop. There is, furthermore, a function call on Line 5 which has to be inlined since it might hide other memory accesses.

Note that, backprop from Rodinia is not a real application but a simplied benchmark that does not interleave calls to malloc and free. Consequently conn[0], conn [START_REF] Ao | Parainsight: An assistant for quantitatively analyzing multi-granularity parallel region[END_REF], . . . , conn[n1] often happen to be laid out contiguous in memory. This is not realistic in practice, though.

By introducing even only a short sequence of calls to malloc and free at the beginning of the benchmark the layout will be non-contiguous as shown in Fig. 2. The reason for this is that malloc gives no guarantees on the placement of allocations. In practice even runtime polyhedral Inria conn [START_REF] Ao | Parainsight: An assistant for quantitatively analyzing multi-granularity parallel region[END_REF] conn [START_REF] Bastoul | Generating loops for scanning polyhedra: Cloog users guide[END_REF] Figure 2: Memory layout for the conn array with n2 = 3 optimizers [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF][START_REF] Simbürger | Polyjit: Polyhedral optimization just in time[END_REF] will only sees random patterns when monitoring memory addresses for the access of conn[k][j] and perform no transformation.

Solution: folding-based analysis

Despite the presence of a non-ane memory access, the above computation kernel presents an interesting opportunity for optimization. Our dynamic analysis detects:

the stride-1 access for conn[k][j] along the outer dimension j;

the absence of dependence along dimension j;

From this information, our back-end suggests to perform a loop interchange, vectorization, and tiling. Applying those transformations lead to a speedup of ×5.3.

The vectorization opportunity is revealed by looking at the scalar evolution [START_REF] Pop | Induction variable analysis with delayed abstractions[END_REF][START_REF] Van Engelen | Ecient symbolic analysis for optimizing compilers[END_REF] of the addresses being accessed, that is, how they change as a function of the values of the iterators k and j. In the case of our example, the addresses used for loading conn [k] as shown in Fig. 2, can be described with the expression 0j + 1k + 12, where 12 is the base address of conn. This is because &conn [k] does not depend on j, and k is incremented by one on each iteration. Note that due to the gap in the layout of conn, the addresses used to access conn[k][j] cannot be described by an ane expression. This is shown in Fig. 3. Our analysis is robust against this irregularity along dimension k and is able to produce the expression 1j + k + 66, where 66 is the base address of the nested array conn[0], and where represents the fact that accesses are not ane along dimension k. However, the obtained expression does indicate that the memory address increases by 1 every iteration of dimension j. We refer to this as a stride-1 access.

The folding algorithm not only discovers the structure of memory accesses, but also the structure of data dependencies in general. In our running example it detects the reduction in sum on line 4 of Fig. 1. It also detects that there is no dependence between the reads on line 4 and the write on line 5. It is worth mentioning that the structure of memory accesses and dependencies are detected separately. The folding algorithm can thus handle cases where accesses are non-ane and dependencies are anes. Here the irregularity of the former does not hinder the folding algorithm from nding the structure of the latter, and vice versa.

The stride-1 access along dimension j allows deducing that SIMD vectorization might be protable. Since j is not the innermost loop it is necessary to perform a loop interchange before vectorizing. That this loop interchange is valid is clear from the absence of dependencies between the two loops. Note that the interchange will require an array expansion of the sum variable along with a new 1-dimensional loop iterating over j to ll the l2 array. Our analysis, like any dynamic approach reasoning on an execution, cannot guarantee that this holds in general, but it can still provide useful feedback. 3 Interface of the Folding-Based Analysis Before describing the core algorithm of the folding-based analysis in Section 4, we introduce its inputs and outputs.

RR

Inputs

The input for the folding-based analysis is provided by our front-end. In order to handle any kind of loops in a uniform way, our front-end inserts canonical iterators in every loop. These iterators start at zero and advance by one every iteration. The front-end is implemented as a plugin for the dynamic binary translator QEMU [START_REF] Bellard | a fast and portable dynamic translator[END_REF][START_REF] Guillon | Program instrumentation with qemu[END_REF]. Even though the front-end analyzes machine code it works at the level of the generic QEMU IR, making it CPU architecture agnostic.

The input of the folding algorithm is composed of streams of two types, one for instructions and one for data dependencies. In the following a static instruction is a machine instruction in the program binary. An instruction instance is one dynamic execution of a static instruction.

A dependence is a pair consisting of an instruction instance that produced a value and another instance consuming it. We call those instances the source and the destination respectively. We note this as source → destination. The folding-based analysis itself does not depend on this and could work at any granularity, such as the source level.

Each input stream has a unique identier Id. An instruction stream is identied by a static instruction, while a stream of data dependencies is identied by a pair of static instructions. The two types of streams have the same overall structure where each entry consists of two elements:

an iteration vector (IV ): a vector made up of the current values of all canonical loop iterators; a label: the denition of the label diers between the two types of streams and is described below.

For a given stream and a given identier, all the IV s span a multi-dimensional space where each entry is a point. Thus, in the following we use the terms entry and point interchangeably. Also, note that IV s arrive in the input stream in lexicographical order.

Instructions An instruction stream for a static instruction Id contains all its instances. The label is a scalar value whose meaning depends on the type of the static instruction. If the instruction is an arithmetic instruction (Cpt) dealing with integers, the label is the integer value representing the result computed by the instruction. If the instruction is a memory access (Mem), the label is the address read or written by the instance. 4 To illustrate the contents of the input stream of instruction instances we again use the example of backprop from Fig. 1. At the binary level, the considered loop-nest contains several instructions that are represented in an abstract C-like fashion in Fig. 4. An excerpt of four instruction streams for this example is shown in Table 1. The IV of each entry is the vector made up of the current values of all canonical loop iterators noted cj and ck in the table.

Dependencies A dependence stream for a pair of static instructions source and destination contains an entry for each pair of instances for these instructions that have a data dependence.

The IV of an entry is the IV of destination, whereas the label is the IV of source. Table 2 shows the dependency input streams of the folding algorithm for the example in Fig. 4. In this example, all the dependencies except I4 → I4 are intra-iteration dependencies.

I1 → I2 I2 → I4 I4 → I4 IV Label IV Label IV Label (cj,ck) (cj',ck') (cj,ck) (cj',ck') (cj,ck) (cj',ck') (0,0) (0,0) (0,0) (0,0) (0,1) (0,1) (0,1) (0,1) (0,

Outputs

The folding algorithm processes each stream independently. For each stream, the nal result of folding is a piecewise linear expression mapping IV s to labels. We refer to this piecewise linear expressions as a label expression. The domain of a label expression contains all the IV s of all points seen in the input stream. We use the terms domain and geometry interchangeably in the following. Each piece of the domain is described by a set of ane inequalities, hence it denes a polyhedron. Section 2.2 already showed two examples of label expressions, e.g., 1j + k + 66.

This is a compact representation of the input stream since it can describe arbitrarily many points in one piece. It also directly exposes regularity in a form that polyhedral optimizers can exploit.

The coecients of a label expression may be either an integer or , as illustrated for the non-ane memory access of backprop shown in Fig. 3. If a coecient is , this indicates that the evolution cannot be expressed as an ane expression along the corresponding dimension.

Instruction For an instruction stream, depending on the type of its corresponding static instruction, the label expression either represents the integer values computed by the instruction or the addresses it accesses. Table 3 illustrates the outputs for the input streams in Table 1 where n2 = 16 and n1 = 42. All instruction instances of every input stream are now described by a single line each. We notice from this table that two of the four instructions have an ane expression where all the coecients are known, that is, they are not . The ane expression of instruction I4 is marked as N/A because it is computing oating point values. Instruction I2 has an ane expression with the coecient for dimension k being , as already discussed.

Nevertheless, the algorithm still outputs the single polyhedron describing the domain for this instruction. It is also worth mentioning that, unlike in this example, the label expression of each instruction can be made up of several pieces. The domain would then be represented as a union of polyhedra.

Dependencies The label expression of a dependency is a piecewise linear expression with multiple outputs. The label expression maps IV s of the consumer instances of the dependence to IV s of the producer instances. That is, given an instruction instance the label expression can be used to determine from which other instruction instances it consumed data. Table 4 illustrates the result of the folding-based analysis for the three dependency input streams in Table 2. All the dependencies of a given input stream are now described by a single line. Each one of these lines states when the dependency between two instruction instances occurs. For example, the last line tells us that the instance (cj, ck) of I4 depends on the instance (cj, ck -1) of itself. As for the output regarding instruction streams, it is worth noting that in this example the domain of all the dependencies is described by a single polyhedron. Nevertheless, in more complex cases these domains can be represented by a union of polyhedra.

Inria Id Polyhedron Label expression (cj,ck) f(cj,ck) I1 → I2 0 ≤ cj ≤ 15, 0 ≤ ck ≤ 42 cj = cj + 0ck, ck = 0cj + ck I2 → I4 0 ≤ cj ≤ 15, 0 ≤ ck ≤ 42 cj = cj + 0ck, ck = 0cj + ck I4 → I4 0 ≤ cj ≤ 15, 1 ≤ ck ≤ 42 cj = cj + 0ck, ck = 0cj + ck -1
Table 4: Output of the folding algorithm for the dependencies stream shown in Table 2 3.3

Using the output

The output of the folding algorithm is intended to be consumed by the back-end of our tool chain leveraging a classic polyhedral optimizer. Such an optimizer requires as input the list of instructions along with their domains and their dependencies. The back-end then searches which re-scheduling transformations can be applied to the instructions under the constraints imposed by the data dependencies.

Before providing dependencies to the back-end, the output stream of dependencies is pruned by removing all the dependencies involving a computation instruction identied as an induction variable. An induction variable is a computation instruction with a label expression where all coecients of all pieces are integers, that is, not . The initial loop iterators are an example of induction variable, e.g., I5 and I6. Removing those instructions serves two purposes. First, induction variables always depend on their value from the previous iteration of the loop they are in.

Consequently their dependencies constrain the execution to be completely sequential. Removing these instructions gives the back-end more freedom and may uncover parallelism or potential for other polyhedral transformations. The second reason for removing induction variables is simply that it reduces the number of instructions the polyhedral back-end has to deal with.

Then, still before providing the dependencies to the optimizer, we must process dependencies having coecients in their label expression. Observe that the fact that some dependencies are not accurately captured by our folding algorithm is not a limitation of the approach, but a choice imposed by polyhedral back-ends which complexity are combinatorial with the size of the polyhedral representation. To that end, we over-approximate those dependencies by imposing a lexicographical ordering over their IV s for the iterators having at least one top coecient. With this order, it is guaranteed that all instances of the producer come before any instances of the consumer that might possibly consume them. For instance, let us assume in our running example that the dependency I4 → I4 is not cj = cj + 0ck, ck = 0cj + ck -1 but cj = cj + 0ck, ck = 0cj + ck: The over-approximated dependency given to the back-end would be cj = cj ∧ ck ≤ ck. Note that, if for the sake of clarity the widening described by our folding algorithm is rather extreme, i.e., the value in our lattice is either an integer or .

Intermediate values such as intervals could be used instead. E.g., see [START_REF] Boulet | Loop parallelization algorithms: From parallelism extraction to code generation[END_REF] for a description of dierent dependency levels.

Finally, the access functions for memory instructions are also given to the polyhedral optimizer so that it can identify opportunities for exposing vectorization and spatial locality. For this it needs information about stride which is given by a non-coecient in the label expression of an instruction accessing memory.

The Folding Algorithm

This section gives an overview of the folding algorithm and then presents the details of the algorithm.

Overview

As stated in the previous section, the folding algorithm processes the stream for each identier separately. For both instruction and dependency streams, the algorithm receives points in a geometrical space as specied by the IV s. In both cases, the main idea of the algorithm is to construct polyhedra from those points. For each polyhedron the algorithm constructs an ane expression describing the label of the points contained in the polyhedron. When receiving the rst point, the algorithm creates a 0-dimensional polyhedron containing only that point. It then tries to grow this polyhedron with the next points, adding dimensions as necessary.

To give an intuition about how the folding algorithm works, let us consider the stream of I2 in Table 1.

Geometric folding

The folding process for I2 is illustrated in Fig. 5. For now we will ignore the construction of the ane expression. As shown, the process leads to the creation of many intermediary polyhedra which are merged as the algorithm executes. The polyhedron P1, a 3 × 3 square, is the nal result of the algorithm. As shown in Fig. 5 the main steps of the algorithm are as follows:

create the 0-dimensional polyhedron P1 when the rst point (cj = 0, ck = 0) is received;

when (cj = 0, ck = 1) is received, P1 absorbs it to become a 1-dimensional polyhedron, that is, a line segment; Inria when (cj = 0, ck = 2) is received, P1 absorbs it;

notice that the loop over ck is completed when point (cj = 1, ck = 0) is received because the iterator of the surrounding loop cj increased. Then create the new 0-dimensional polyhedron P4;

P4 absorbs (cj = 1, ck = 1) to become a 1-dimensional polyhedron and then absorbs (cj = 1, ck = 2) (not shown in Fig. 5);

notice that the loop over ck is completed when point (cj = 2, ck = 0) is received. P1 absorbs P4 along dimension cj. Then create the new 0-dimensional polyhedron P7;

P7 absorbs (cj = 2, ck = 1) to become a 1-dimensional polyhedron and then absorbs (cj = 2, ck = 2) (not shown in Fig. 5);

P1 absorbs P7 and becomes the nal 3 × 3 square.

The geometric folding works exactly the same for dependencies as illustrated above for instructions. The only dierence is the semantic of the reconstructed union of polyhedra. In the case of an instruction, this union denes when the instruction is executed. For a dependency it tells when the dependency occurs from the point of view of the destination.

Label folding

In the previous section we ignored the folding of the labels associated with each point in the input stream. Nevertheless, this label folding takes place at the same time as geometric folding.

It is also performed in a streaming fashion. In the context of label folding, the symbol ⊥ denotes a coecient that has not yet been determined because the loop has not yet iterated along the dimension associated with that coecient. As shown in Fig. 5 the label folding proceeds as follows:

create f 1(cj, ck) = ⊥cj + ⊥ck + 67 when point (cj = 0, ck = 0) with label 67 is received;

update f 1 to ⊥cj + 4ck + 67 when P1 absorbs (cj = 0, ck = 1) with label 71 is received because ck advanced by 1 and 71 -67 = 4;

check if f 1(cj, ck) = ⊥cj + 4ck + 67 is valid when P1 absorbs (cj = 0, ck = 2) with label 77. It is not the case, so update f 1 to ⊥cj + ck + 67; repeat the steps above for P4 and get f 4(cj, ck) = ⊥cj + ck + 68 (not shown in Fig. 5); update f 1 to f 1(cj, ck) = 1cj + ck + 67 when P1 absorbs P4 because cj advanced by 1 and 68 -67 = 1;

check whether f 1(cj, ck) = 1cj + ck + 67 is compatible with f 7(cj, ck) = ⊥cj + ck + 69, when P7 absorbs P1 to get the nal 3 × 3 square. It is the case.

The algorithm that folds the labels of a dependency is the same as the one described above for the label of an instruction. It is just applied individually for each scalar value in the label vector, that is, each component of the IV of the source of the dependency.

The algorithm

This section introduces the structure of the main algorithm itself and then explains its subcomponents.

RR n°9244 4.2.1 Main folding function

The main function is shown in Algorithm 1. As explained in Section 3, this main function, is applied to each input stream separately. To handle real-life applications, where input streams are huge, the algorithm works in a streaming fashion (Line 10). It is not necessary to have the whole input available at once. The output is also emitted as a stream. The main principle of the algorithm, as depicted in the example in Fig. 5, consists of maintaining a worklist of intermediate polyhedra per dimension. The intermediate polyhedra then grow by absorbing other polyhedra.

Note that a d-dimensional polyhedron can only absorb (d -1)-dimensional polyhedra.

Representation for polyhedra The folding algorithm only produces bounded convex polyhedra (polytope). Internally, these polyhedra are represented by their extreme points (vertices).

Also, all the polyhedra have edges whose slopes are always in {-1, 0, 1}. We call those, elementary polyhedra, and dene them formally using the following recursive denition:

an elementary 0-dimensional polyhedron is a polyhedron made of a single point;

an elementary d-dimensional polyhedron is a convex polyhedron with 2 d extreme points such that: 1. All its extreme points must have identical coordinates in dimensions higher than d; 2. Its lower and upper faces must themselves be (d -1)-elementary polyhedra; 3. The edges connecting the lower and upper faces can be expressed as k S. Where k ∈ N * and S, the slope vector of the edge, is a vector where all components are either -1, 0, or +1.

A polyhedron is degenerate on a given dimension if all its extreme points have the same coordinate for that dimension, that is, it has zero width in that dimension. The elementary polyhedra produced by the folding algorithm may be degenerate on one or more dimensions.

Producing only elementary polyhedra allows the absorption process described below to work is the extreme point from which the absorption is done. This point is the lexicographically rst point of the upper face of abso.

Analysis steps When a point is received, the algorithm rst processes the innermost dimension (numbered 1). Then, for each loop (but for the outermost one) that completes in the instrumented code, the algorithm processes its enclosing dimension. In other words, if the innermost loop nishes, the algorithm processes dimension d = 2; if its enclosing loop nishes, it processes dimension d = 3; etc. In Line 17, process_dims represents that set of dimensions to be processed.

Inria

Before processing the dierent dimensions, the current point is added into absorbers[0]

(Line 14). This state is only transient, because as soon as the innermost dimension is processed, the point will be promoted into eps_2_to_be_absorbed [START_REF] Ao | Parainsight: An assistant for quantitatively analyzing multi-granularity parallel region[END_REF] (Line 21). Then, for each dimension d of process_dims (processed from inner to outer), three steps are performed. The rst step (Lines 20 to 21) promotes all polyhedra in absorbers[d-1] into eps_2_to-_be_absorbed [d]. Because dimensions are processed in increasing order, that is, from innermost to outermost, when processing dimension d we are sure that absorbers[d-1] have already absorbed all the (d-2)-dimensional polyhedra it could. This promotion to d-dimensional degenerate polyhedra allows them to be absorbed in the next step by the d-dimensional polyhedra already in absorbers [d].

In the second step (Lines 24 to 39), polyhedra from absorbers[d] try to absorb polyhedra in eps_2_to_be_absorbed [d]. For absorption to be possible, the polyhedra should be geometrically compatible (Line 30) and their label expressions should match (Line 31) as described in Section 4.2.1 and Section 4.2.2. If a polyhedron in absorbers [d] does not absorb any other polyhedron, then it will never grow again along dimension d. As a consequence, it is promoted into eps_2_to_be_absorbed[d+1] (Line 39). This promotion also transforms the d-dimensional polyhedron into a (d + 1)-dimensional degenerate polyhedron.

The third and last step (Lines 43 to 44) promotes all the d-dimensional polyhedra in eps_2_to-_be_absorbed[d] that have not been absorbed. Since those polyhedra will never be absorbed again in dimension d, they are moved to the absorbers[d] list so that they will have a chance to themselves absorb other polyhedra next time dimension d is processed.

During the execution of the algorithm, a polyhedron is retired when it is promoted to the dimension above the maximum dimension of the space, e.g., 3 for an instruction in a 2D loop nest. When the stream is nished, all remaining non-retired polyhedra are also retired. Retired polyhedra are written to the output stream and do not consume memory anymore. This is safe since we know that they will never grow anymore.

Absorption

As stated in Section 4.1, the second step of the folding algorithm grows polyhedra by letting them absorb each other. A d-dimensional polyhedron searches for candidates to absorb by checking if the rst point of its upper face, called the corner (Line 27), touches the anchor of any other (d -1)-dimensional polyhedron. This search is done by adding the search vectors v to the coordinates of the corner and performing a lookup in eps_2_to_be_absorbed [d] to see if there is a polyhedron at this position (Line 28). Once a candidate has been found, the algorithm must check that all other extreme points also match (Line 30). Which search vectors are used for this lookup and how geometric matching is checked depends on whether the absorber is degenerate in d or not. If the absorber is degenerate we call this a polyhedra merge. An example of this is when P1 absorbs P4 in Fig. 5. The second case, a polyhedra extension, occurs when the absorber is not degenerate, as seen for example when P1 absorbs P7.

Polyhedra merge In this case, the d-dimensional absorber polyhedron is degenerate on dimension d. Hence, it has no slopes yet on that dimension. As a consequence the set of search vectors used to nd candidates are all possible slope vectors where the value of the i t h component is a)

0 if i > d; b) 1 if i = d; c) in {0, -1, +1} if i < d.
A polyhedra merge is legal if every extreme point P of the upper face of the absorber can be connected to the corresponding extreme point of the polyhedron to be absorbed using any slope vector strictly greater than P . Note that there are always exactly 3 d-1 such slope vectors and that they can be pre-computed. After absorption the resulting polyhedron will no longer be degenerate in d. Its lower face will be the lower face RR n°9244 of the original absorbing polyhedron and the upper face of the absorbing polyhedron will be the extreme points of the (d -1)-dimensional absorbed polyhedron. The slopes of edges between the faces of the new polyhedron will be exactly those search vectors used to connect the extreme points.

Polyhedra extension In this case the absorber is a non-degenerate d-dimensional polyhedron.

Hence, the absorber already has slopes for all its edges. When looking for candidates to absorb, there is only one search vector, the slope vector connecting the oldest extreme point of the lower face to that of the upper face. To check if it is legal for the absorber to absorb the candidate it suces to verify whether the extreme points of the two polyhedra can be connected using the slopes of the absorber.

After an absorber polyhedron has found a candidate with extreme points that match, the algorithm checks whether the labels also match. This is done by the has_compat_label function described in the next section. If both the geometry and the label are compatible, then the absorption is performed by a call to absorbs.

Compatibility and update of label expressions

The data structure used for label expressions is shown in Fig. 6. num_dimensions is the number of loops enclosing the static instruction or the destination instruction associated with the input stream. These coordinates are used when coecients are updated. Note that once a coecient has been updated from an unknown to a known value, it is never updated again except to be set to .

Label_Function: int num_dimensions int[num_dimensions + 1] init_point int[num_dimensions + 1] coeffs coeff_t[num_dimensions + 1] coeff_types

At creation time, coeff[0] is given the value associated with the initial point. As long as there are some ⊥ coecients, coeff[0] contains the remaining amount contributed by unknown coecients. We refer to coeff[0] as the remaining value in the following. This remaining value is updated whenever a coecient is updated. When all coecients are known, the remaining value represents the constant coecient of the ane expression.

The two polyhedra involved in a compatibility check along dimension d may be degenerate on one or more dimensions, including the d th one. As a consequence, the check may be faced with ane expressions where some coecients are ⊥. In the following, we note the label expression of the absorbing polyhedron as f_abs, and that of the polyhedron to be absorbed as f_to_be_abs.

Inria

We notice that the polyhedron to be absorbed is always degenerate on dimension d, as stated in Section 4.2.1. Hence, f_to_be_abs.coeff_types[d] = ⊥.

All dimensions below are known For illustrative purposes we rst cover the simplied case where all dimensions below d are known for the two label expressions. Here the compatibility check, is straightforward. First the function has_compat_label veries that all coecients for dimensions from 1 to d -1 are the same. If this is not the case the two label expressions are incompatible, and it returns false.

Otherwise, the check may be faced with two cases corresponding to the two dierent absorption cases described in Section 4.2.2. In the polyhedra merge case, where the absorber polyhedron is degenerate on dimension d, that is, f_abs.coeff_types[d] = ⊥, the check always succeeds and the has_compat_label function returns true. Indeed, by setting the proper coecient for dimension d and by updating the remaining value, it is always possible to make the two expressions compatible as shown by the update_label_dims_known function in Algorithm 2. The new coecient is equal to the dierence of remaining values (Line 6). Note that, in general we would also have to divide the new coecient by the progress made along dimension d. However, because absorption guarantees that the two polyhedra whose label expressions are being merged touch each other the progress is always equal to 1. Finally, the remaining value is decreased by the eective contribution of the new coecient taking into account the d th coordinate of the initial point (Line 10).

In the polyhedra extension case, the absorber polyhedron is not degenerate on dimension d. Its ane expression already has a value computed for the coecient on dimension d Then f_abs.coeff_types[d] = ⊥ and nothing needs to be updated. The compatibility check must only ensure that this coecient is compatible with f_to_be_abs. Algorithm 2 shows the has_compat_label_dims_known function implementing this check. First, it computes the contribution of the known coecient of f_abs into f_to_be_abs using the initial point of f_to_be_abs (Line 14). Then, the check subtracts this contribution from the remaining value of f_to_be_abs to compute its new remaining value. For the check to return true, this new remaining value must be equal to the remaining value of f_abs (Line 16).

General case In the general case the two polyhedra may be degenerate for some dimensions below d. This happens if a dimension only iterates once. The compatibility check described above must take this into account.

Function has_compat_label_general in Algorithm 3 shows the general compatibility check. If f_abs does not have a coecient set for dimension d, the check always succeeds as in the particular case previously described. The check works by comparing the coecients of both expressions for all the dimensions from 1 to d. If both coecients for a dimension are known they must be the same or the check fails (Line 9). If one is known and not the other (Line 12 and Line 14), then the function increments the total contribution coming from the other expression for the expression having the unknown coecient. At the end of the loop, the check ensures that the coecient for dimension d in f_abs is compatible with f_to_be_abs. This check relies on the total contribution variables incremented during the loop to ensure that the two expressions still produce the same value after merging.

In case they are compatible, the new coecients, that is, the one on dimension d and potentially others, and the new remaining value for the expression of the absorber are computed by the same principles as the ones performed by the has_compat_label_dims_known.

Label widening As shown by the backprop example, the folding algorithm must be capable of identifying labels that are ane on some dimensions and not on others. To that end, the RR n°9244 algorithm has a mechanism called label widening enabling it to skip the matching of labels on a per dimension basis. If the compatibility check between two coecients fails, then instead of returning false (Line 11 in Algorithm 3), the coecient is set to and true is returned instead. The absorption can still happen, even if the labels of the two polyhedra are not fully compatible. The resulting polyhedron is no longer a fully accurate representation of the input stream. Nevertheless, this mechanism allows the folding algorithm to handle real life applications without a perfect ane behavior. The name label widening stems from the fact that in the case of dependencies it widens the label expressions from strict equalities to inequalities, as shown in Section 3.3.

The integration of this feature into Algorithm 2 and Algorithm 3 is straightforward. A coecient is compatible with any other coecient, and when performing absorption, any such coecient in one of the two label expressions leads to a coecient in the updated expression.

The label widening mechanism is crucial for the label expressions of instructions because is a clear indicator that a memory accesses is not ane along a dimension. For dependencies it simply reduces the size of the output given to the back-end by reducing the number of produced pieces.

Geometric give up

Even with the label widening mechanism described above, some applications may lead to the creation of a huge number of polyhedra. This happens when the geometry of instructions and dependencies are not ane. In the worst case, the folding algorithm creates one polyhedron for each dynamic instruction and for each dynamic dependency.

To mitigate this issue, the folding algorithm has another global option called geometric giveup. This options allows dening an upper limit on the number of intermediate polyhedra. Remember that an intermediate polyhedron is a polyhedron in one of the worklists that can still grow by absorbing other polyhedra. Before creating a new polyhedron (Line 14), the algorithm checks if the number of intermediate polyhedra exceeds the threshold. If so, then the associated input stream is marked as give up. Once a stream has been marked as give up, nearly all information heretofore collected for it is discarded. First, all intermediate polyhedra and also the polyhedra that have already been retired are discarded. Furthermore, all coecients for all outputs of the label expression are set to , that is, a geometric give up implies giving up on all dimensions of the label expression. The only information that is retained for the stream is the maximum coordinates seen in the IV s of any point. From then on every time a new point is received for the given up stream, the folding algorithm previously described is skipped. Instead, only the maximum coordinates seen are update as necessary for every points.

The nal geometry emitted for a give up stream is simply an hyperrectangle that starts at the origin and extends to the maximum coordinates seen in the IV s of any point of the input.

In other words, the geometry of the input stream is over-approximated by a large polyhedron.

Experimental Results

This section applies our analysis to a full benchmark suite to demonstrate the scalability of the folding algorithm and show that it extracts rich information for optimization.

Experimental setup We use the latest revision, 3.1, of the Rodinia benchmark suite [START_REF] Che | Rodinia: A benchmark suite for heterogeneous computing[END_REF][START_REF] Che | A characterization of the rodinia benchmark suite with comparison to contemporary cmp workloads[END_REF].

All measurements and experiments where performed on a Xeon Ivy Bridge CPU with two 6 core CPUs, each running at 2.1GHz. As the front-end producing the IV s and labels does not support multithreaded applications yet, each benchmark is run with a single thread. All Inria benchmarks where compiled using GCC 8.1.1. Since QEMU, which the front-end is based on, currently cannot handle newer AVX instructions we used the compiler ags -g -O2 -msse3. For the speedup measurements of backprop mentioned in Section 2.2 we used the Intel icc 18.0.3 compiler and the ags -Ofast -march=native -mtune=native.

Note that the instructions in our experiments are real X86 machine instructions. Many X86 instructions both read or write memory and perform computations at the same time. As a consequence the instructions streams that form the input of the folding algorithm are actually more complicated than the ones presented in Section 3.1 and Table 1 in a simplied way for clarity purposes. In reality the label of an instruction can have multiple values to account both for the addresses accessed and the values produced. The label expressions for instructions thus potentially have multiple outputs as well, just like those for dependencies.

Table 5 gives statistics on the size and precision of the output of four versions of the folding algorithm. F is the basic algorithm as described in Section 4, with label widening for instructions and without for dependencies. F W is the algorithm with label widening for both instructions and dependencies. F GG is the same as F but with geometric give up. F GG,W is the same as F W but with geometric give up. The treshold for the geometric give up was set to allow 4d + 1 intermediate polyhedra in each d dimensional space. That is, enough for the ane expression constructed to be made up of up to four d dimensional pieces.

For each algorithm we report the following statistics. #P is the number of polyhedra in the output stream. For dependencies, %A is the number of dependence instances that where in an ane piece of the label expression. A piece of the label expression is considered ane if it has no coecient. This column is omitted for algorithm F since by construction it always contains 100%. Similarly for instructions, %A is the number of instruction instances that where in an ane piece of the label expression. A piece of the label expression of a static instruction is considered ane if it either: (a) does not perform a memory access, or (b) has no coecient in its memory access function. #MP l is the maximum number of intermediate polyhedra live at any moment of the execution, indicating the memory usage of the algorithm.

The remaining columns in the table are as follows. Input Size shows the total number of entries in all dependency and instruction input streams. Optim shows a very brief outline of the optimization feedback given by our polyhedral back-end using the output of F GG,W [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF]. In this column T nD indicates that the back-end has found that n dimensional tiling was possible. P indicates that the back-end has detected parallelism that can be exploited using threads.

V indicates that the back-end has detected potential for vectorization. Note that the entire feedback of the tool is immensely richer and more elaborated [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF], this column gives only a simplied summary.

Finally, note the numbers reported in Table 5 corresponds to applying the folding-based analysis on the hot region of each benchmark, we have ltered out the phases where the benchmarks read their input or write their output. This hot region often involves numerous function calls [START_REF] Gruber | Data-ow/dependence proling for structured transformations[END_REF].

Discussion of the results Since the polyhedral optimization performed in the back-end is an exponential problem it is crucial that the output of the folding-based analysis is of tractable size. Table 5 clearly shows that F GG and F GG,W produce drastically smaller outputs than the other two versions. As indicated by the %A column, F GG,W is roughly as precise as F GG , but produces an even smaller output. In fact only the output of F GG,W is small enough for the back-end to handle.

Since Rodinia is a benchmark suite designed to exploit multi-core parallelism each benchmark contains at least one parallel loop. As seen in column Optim the folding-based analysis clearly detects this parallelism across the entire suite, even in the presence of may-alias dependencies in the source code. We also nd that there is tiling potential across Rodinia. Note that streamcluster, the least ane of all benchmarks, exhausted memory in the polyhedral back-end and therefore no result is displayed. Benchmark mummergpu is not included in the results since it contains CUDA code and the front-end can only instrument code run on the CPU.

RR n°9244

Related Work

Integer linear algebra is a natural formalism for representing the computation space of a loop nest. The polyhedral framework [START_REF] Feautrier | Polyhedron model[END_REF] leverages, among others, operators on polyhedrons, enumeration for code generation [START_REF] Bastoul | Generating loops for scanning polyhedra: Cloog users guide[END_REF], and parametric integer linear programming [START_REF] Feautrier | Parametric integer programming[END_REF] for dependence analysis [START_REF] Collard | Fuzzy array dataow analysis[END_REF]. Historically, it has been designed to work on restricted programming languages, and was used as a framework to perform source-to-source transformations. More recently, eorts have been made to integrate the technology in mainstream compilers with GCC-Graphite [START_REF] Trifunovic | GRAPHITE Two Years After: First Lessons Learned From Real-World Polyhedral Compilation[END_REF] LLVM-Polly [START_REF] Grosser | Polly -performing polyhedral pptimizations on a low-level intermediate representation[END_REF]. The set of loop transformations that the polyhedral model can perform is wide and covers most of the important ones for exposing locality and parallelism to improve performance [START_REF] Bondhugula | A practical automatic polyhedral program optimization system[END_REF].

Dynamic data ow/dependence analysis is a technique typically used to provide feedback to the programmer, e.g., about the existence or absence of dependences along loops. The detection of parallelism along canonical directions, such as vectorization, has been particularly investigated [START_REF] Li | Discopop: A proling tool to identify parallelization opportunities[END_REF][START_REF] Butt | Runtime analysis of application binaries for function level parallelism potential using qemu[END_REF][START_REF] Wang | Integrating proledriven parallelism detection and machine-learning-based mapping[END_REF][START_REF] Ketterlin | Proling data-dependence to assist parallelization: Framework, scope, and optimization[END_REF][START_REF] Ao | Parainsight: An assistant for quantitatively analyzing multi-granularity parallel region[END_REF][START_REF] Kim | Prospector: A dynamic data-dependence proler to help parallel programming[END_REF][START_REF] Faxén | Embla -data dependence proling for parallel programming[END_REF][START_REF] Vandierendonck | The paralax infrastructure: automatic parallelization with a helping hand[END_REF][START_REF] Tournavitis | Semi-automatic extraction and exploitation of hierarchical pipeline parallelism using proling information[END_REF], as it requires only relatively localized information.

Another use case is the evaluation of eective reuse [START_REF] Marin | Miami: A framework for application performance diagnosis[END_REF][START_REF] Liu | Pinpointing data locality problems using datacentric analysis[END_REF][START_REF] Beyls | Discovery of locality-improving refactorings by reuse path analysis[END_REF][START_REF] Berg | Fast data-locality proling of native execution[END_REF] with the objective of pinpointing data-locality problems. Like us, with the objective of gathering a more global dependence information, Redux [START_REF] Nethercote | Redux: A dynamic dataow tracer[END_REF] builds a complete extended dynamic dependence graph from binary level programs. The paper concludes with a negative result. Because of its inability to compress the produced graph it is only able to handle very small non-realistic programs.

Among the existing trace compression algorithms, two specialize in extracting a polyhedral representation from input streams [START_REF] Ketterlin | Prediction and trace compression of data access addresses through nested loop recognition[END_REF][START_REF] Rodríguez | Trace-based ane reconstruction of codes[END_REF]. However, although they excel in rebuilding a polyhedral representation for a purely ane trace, they suer inherent limitations for (even partially) non-ane traces. They share the idea of using pattern matching with ane expressions with our folding algorithm but do not exploit the geometric information provided by the IV s. Nongeometry-based approaches require a nite window of points under consideration. Unfortunately, this forces a trade-o between speed for and quality of the output when choosing the size of this window. For perfectly regular programs a small window can be used, making the algorithms Inria very ecient. In that simple case using a geometric approach does not make much dierence.

With d the dimension of the iteration space and n the number of points, the complexity of both non-geometric approaches is O(2 d n). However, in the context of proling large non-fully ane programs, none of these two existing approaches can be used. The complexity of the nested loop recognition algorithm of Ketterlin et al. [START_REF] Ketterlin | Prediction and trace compression of data access addresses through nested loop recognition[END_REF] increases quadratically with a parameter k that bounds the size of the window. If k is smaller than the amount of irregularity along the innermost dimension, it is not able to capture the regularity, and thus compress, along outer dimensions. Hence, to be as ecient as our folding algorithm on our backprop example, for most execution instances, k would have to be bigger than 10 4 . The complexity of the ane recognition algorithm of Rodriguez et al. [START_REF] Rodríguez | Trace-based ane reconstruction of codes[END_REF] increases exponentially with the number of irregularities. So in practice, it has to give up even for nearly ane traces.

Similarly to us, existing runtime polyhedral optimizers [START_REF] Simbürger | Polyjit: Polyhedral optimization just in time[END_REF][START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF] use runtime information to create a polyhedral representation of a program. PolyJIT [START_REF] Simbürger | Polyjit: Polyhedral optimization just in time[END_REF] focuses on handling programs that do not t the polyhedral model statically because of memory accesses, loop bounds and conditionals that are described by quadratic functions involving parameters. Apollo [START_REF] Martinez Caamaño | Full runtime polyhedral optimizing loop transformations with the generation, instantiation, and scheduling of code-bones[END_REF] handle this case and many others preventing static polyhedral optimizers from operating. Compared to our analysis, PolyJIT focuses on identifying 100% ane programs which may be rare in practice for many reasons such as the memory allocation concerns pointed out for backprop.

Apollo proposes a tube mechanisms [START_REF] Sukumaran-Rajam | The polyhedral model of nonlinear loops[END_REF] which allows the handling of programs with quasi-ane memory accesses. Even with this last extension, on the illustrative example of backprop, Apollo will, as opposed to our analysis, neither manage to over-approximate the non constant stride along the inner-most dimension as soon as the stride distance is greater than a given threshold, nor detect the stride of 1 along the outermost dimension. Also, it is worth mentioning that, as for the example of backprop, a program might show ane dependencies while having nonane memory accesses. Contrary to our analysis front-end, which tracks both separately, Apollo only traces memory accesses and then recomputes the dependencies from them. Consequently, Apollo has to give up completely here while we can detect an accurate polyhedral representation of dependencies.

Conclusion and Perspectives

We have presented a folding algorithm able to create a polyhedral representation of a program from its execution trace. Based on a geometric approach, our algorithm scales to real-life applications by safely over-approximating the dependencies that do not t the polyhedral model while still recovering precise information for those that do. From what we observed on the Rodinia benchmark, large regions of programs that are perfectly ane and suitable as-is for classic polyhedral optimizers are rare. Thanks to our over-approximation mechanisms, we are nevertheless still able to create a polyhedral representation for these programs that can be given to a classic polyhedral optimizer.

Regarding the perspectives opened by this work, we are already working in two directions that will allow handling more programs. The rst one consists of adding new dimensions not present in the program to our representation. Said dierently, an instruction contained in a 2-dimensional loop nest in the program could be represented by a 3-dimensional polyhedron. This mechanism, already at work in trace compression algorithms [START_REF] Ketterlin | Prediction and trace compression of data access addresses through nested loop recognition[END_REF][START_REF] Rodríguez | Trace-based ane reconstruction of codes[END_REF] will allow our analysis to handle tiled stencil computations and programs where 2-dimensional arrays are traversed by linearized 1dimensional loops. The second extension we want to investigate is a clever mechanism for the activation of widening for dependency label expressions. We are planning to replace the existing user controlled global option with an adaptive mechanism that automatically activates widening RR n°9244 # Per dimension list of absorber polyhedra. <int, poly_list_t> absorbers # Per dimension dictionary mapping # extreme points to polyhedra to be absorbed <int, <point_t, poly_t>> eps_2_to_be_absorbed # While we have points while

Figure 4 :

 4 Figure 4: C-like binary version for the code of Fig. 1

Figure 5 :

 5 Figure 5: Folding process for the input stream in Fig. 3 considering only three points in both dimensions.

  in no more than O(3 d-1 ) time. The choice of producing only such polyhedra is also motivated by the nature of the input streams that we want to process. The front-end we use to feed the folding algorithm always produces IV s starting at zero and only ever advancing by one. Hence, elementary polyhedra are able to represent the behaviour of most regular loops.Data structures Folding works on spaces with a xed number of dimensions d, that is, the dimensionality of the corresponding IV s. The state of the folding algorithm is contained in two dictionaries. The rst one, absorbers (Line 2), contains a list of intermediate polyhedra for each dimension. absorbers[d] only contains d-dimensional, potentially degenerate, polyhedra. The polyhedra in absorbers[d] are those that can still grow along dimension d by absorbing (d -1)-dimensional polyhedra. Those (d -1)-dimensional polyhedra are stored in eps_2_to_-be_absorbed[d] (Line 6). The keys of the dictionary eps_2_to_be_absorbed[d] are the lexicographically rst, extreme points of the polyhedra to be absorbed. This point, which we name the anchor, is used to uniquely identify the absorbed polyhedron. The abso.upper_left(Line 27)

Figure 6 :

 6 Figure 6: The data structure used to represent label expressions

Algorithm 2 :

 2 (True): pt = wait_next_point() if point == end_of_stream: break # Put current point in absorbers[0] absorbers[0].insert(new Polyhedron(point)) # for each dimension d such that d=1 or loop d-1 completed for d in process_dims(pt): # Step 1: promote absorbers[d-1] -> eps_2_to_be_absorbed[d] for p in absorbers[d-1]: p.move(absorbers[d-1], eps_2_to_be_absorbed[d]) # Step 2: absorbers[d] try to absorb eps_2_to_be_absorbed[d] for abso in absorbers[d]: absorbed = False for v in abso.search_vectors: corner = abso.upper_left to_be_abs = eps_2_to_be_absorbed[d][corner + v] if to_be_abs != None: if (abso.can_absorb(to_be_abs, d) and abso.has_compat_label(to_be_abs, d)): abso.absorbs(to_be_abs, d) absorbed = True break if not absorbed: # abso will never absorb anyone along d, # then promote it in the next dimension abso.move(absorbers[d], eps_2_to_be_absorbed[d+1]) # Step 3: promote all of remaining # eps_2_to_be_absorbed[d] -> absorbers[d] for not_abs in eps_2_to_be_absorbed[d].values: not_abs.move(eps_to_2_to_be_absorbed[d], absorbers[d]) # Stream finished, flush all pending polyhedra flush_pending_polyhedra() Algorithm 1: The main folding algorithm RR n°9244 # Update coefficient for dimension d and remaining # value of f_abs. No need to update f_to_be_abs # because it will be thrown after absorption def update_label_dims_known(f_abs, f_to_be_abs, d): # Update of coefficient new_coeff = f_to_be_abs.coeffs[0] -f_abs.coeffs[0] f_abs.coeffs[d] = new_coeff # Update of remaining value new_coeff_contrib = new_coeff * f_abs.init_point[d] f_abs.coeffs[0] = f_abs.coeffs[0] -new_coeff_contrib # Check in case absorber already has a coeff for d def has_compat_label_dims_known(f_abs, f_to_be_abs, d): new_coeff_contrib = f_abs.coeffs[d] * f_to_be_abs.coeffs[d] new_remain = f_to_be_abs.coeffs[0] -new_coeff_contrib return new_remain == f_abs.coeffs[0] Update and compatibility check when all dimensions below d are known def has_compat_label_general(f_abs, f_to_be_abs, d): if f_abs.coeff_types[d] == ⊥: return True abs_diff = 0 to_be_abs_diff = 0 for q in [1, d]: abs_t = f_abs.coeff_types[q] to_be_abs_t = f_to_be_abs.coeff_types[q] if abs_t != ⊥ and to_be_abs_t != ⊥: if f_abs.coeffs[q] != f_to_be_abs.coeffs[q]: return False if abs_t == ⊥ and to_be_abs_t != ⊥: abs_diff += f_to_be_abs.coeffs[q] * f_abs.init_point[d] if abs_t != ⊥ and to_be_abs_t == ⊥: to_be_abs_diff += f_abs.coeffs[q]*f_to_be_abs.init_point[d] return f_abs.coeff[0] -abs_diff == f_to_be_abs.coeff[0] -to_be_abs_diff Algorithm 3: General compatibility check Inria RESEARCH CENTRE GRENOBLE -RHÔNE-ALPES Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Table 1 :

 1 Instruction input streams from example in Fig.

	. .	. . .	. . .	. . .	. . .	. . .	. . .	. . .

Table 3 :

 3 Output of the folding algorithm for the instructions stream shown in Table1with n2 = 16 and n1 = 42

Table 5 :

 5 Size #P #MP l #P %A #MP l #P %A #MP l #P %A #MP l Size #P %A #MP l #P %A #MP l Evaluation of the folding algorithm

	Benchmark	Input	F	Dependencies F W	F GG	F GG,W Input	Instructions F W	F GG,W	Optim
	backprop bfs b+tree cfd heartwall hotspot hotspot3D kmeans lavaMD leukocyte lud myocyte nn nw particlefilter pathfinder srad_v1 srad_v2 streamcluster	19M 160 5M 903K 95M 91K 782M 530 33G 3K 19M 11K 235M 168 1G 135 1G 7K 5G 516K 89M 2K 4M 5K 782K 124 217M 301 3G 5K 74M 35 3G 250 1G 276 2G 1M		385 160 100% 385 160 100% 385 160 100% 385 15M 140 99% 304 140 99% 304 T 2D, P, V 965K 874K 93% 951K 74 31% 772 70 31% 772 4M 520K 82% 472K 38 51% 367 T 2D, P 390K 86K 99% 336K 113 99% 3K 113 99% 3K 61M 50K 90% 153K 160 89% 1K T 3D, P, V 1K 525 98% 1K 530 100% 1K 525 98% 1K 498M 332 100% 961 332 100% 961 T 3D, P, V 8K 2K 90% 6K 1K 10% 5K 1K 10% 5K 18G 1K 69% 3K 1K 9% 3K T 5D, P 22K 10K 95% 21K 785 0% 6K 785 0% 6K 11M 6K 71% 13K 520 0% 3K T 2D, P 1K 162 91% 1K 168 100% 1K 162 91% 1K 183M 84 85% 782 84 85% 782 T 3D, P 477 131 99% 472 135 100% 477 131 99% 472 911M 82 95% 281 82 95% 281 T 4D, P, V 2K 7K 94% 2K 7K 100% 2K 7K 94% 2K 923M 4K 71% 1K 4K 71% 1K T 3D, P 161K 514K 99% 113K 162 99% 66K 162 99% 65K 2G 355K 84% 72K 128 84% 40K T 3D, P, V 1K 2K 98% 1K 2K 98% 1K 2K 98% 1K 51M 1K 97% 864 1K 97% 864 T 3D, P 9K 5K 100% 9K 5K 100% 9K 5K 100% 9K 3M 3K 99% 4K 3K 99% 4K T 1D, P, V 242 124 100% 211 124 100% 241 124 100% 211 855K 160 100% 189 160 100% 189 T 1D, P 1K 296 99% 1K 301 100% 1K 296 99% 1K 111M 155 100% 555 155 100% 555 T 2D, P, V 92K 3K 99% 2K 550 8% 2K 541 8% 2K 2G 2K 99% 1K 474 11% 1K T 2D, P, V 139 35 100% 135 35 100% 139 35 100% 135 42M 24 61% 116 24 61% 116 T 2D, P 851 242 94% 824 250 100% 851 242 94% 824 2G 179 93% 531 179 93% 531 T 2D, P 811 268 97% 791 276 100% 811 268 97% 791 721M 204 93% 493 204 93% 493 T 2D, P 1M 1M 85% 1M 8K 85% 13K 6K 85% 12K 1G 611K 71% 618K 3K 71% 6K -

as needed. For example, the option could be activated when the number of polyhedra used to represent a given instruction or dependency is becoming too large. This would allow having a trade-o between the accuracy and the size of the output of the folding algorithm.


Razanne Issa

David James Sherman

Deriving Gene Associations for Metabolic Models

Keywords: metabolic networks, comparative genomics, inference réseaux métaboliques, génomique comparée, inférence {G1}

We define a formal procedure for inferring gene-protein-reaction (GPR) relations from complete metabolic models, using a logical representation of knowledge and a small set of inference rules. We show that different use cases of metabolic models requires difference GPR relations. Three examples from the yeast Saccharomyces cerevisae illustrate the procedure and demonstrate its usefulness.

M 0 M 1 m 0 , a {G0} 0 m 1 , a {G1} 1 m 1 , a {G1} 1 (2 )
(3 )

G 0 h • • • G 1 ∼ = (1 )
Figure 1: Incremental modeling: we (1 ) define a tuple m, a {G0} that annotates a minimal enzymatic model m with annotations a 0 defined over a set of variables G, and (2 ) project the tuple from a complete model M . The objective is that (3 ) the inferred model m 1 , a 1 Introduction Gene Associations or Gene-Protein-Reaction relations (GPR) [RVSP03, NvEF + 06] in metabolic models record knowledge about the direct and indirect relations between genes and biochemical reactions. Whether written as boolean formulas or as collections of gene, protein, reaction) triplets, gene associations define boolean functions that assert whether a reaction may take place. One widespread use of gene associations is in metabolic model inference for new target organisms. For such organisms there is typically a genome but very little systematic experimental data. The genes in the genome are used to evaluate the gene associations for every reaction in a reaction database; a reaction is included in the inferred model it its gene association is satisfied.

Gene associations play an essential role in incremental model inference (figure 1) as defined in [START_REF] Issa | symbolic analysis and inference of metabolic models[END_REF]. The incremental method is designed for cases when the target organism is closely related to an organism for which a minimal model m, a {G0} , derived (or derivable) from a complete model M 0 , already exists. The gene associations a {G0} annotate the metabolic model m, and are essential in inferring the target model m 1 , a {G1} 1 . In the general case, however, gene associations are an annotation determined by a human curator on the basis of existing knowledge about the biochemistry and gene regulation of the target organism. While expert creation is the of course the ultimate criterion for correctness, it would be useful to also have a formal criterion defined in terms of the knowledge base.

In this report we describe how a generalized form of gene associations can be used to annotate a minimal enzymatic model, in order to infer specific metabolic models for related organisms based on genomic comparisons. The gene associations will implicitly represent knowledge derived from a hypothetical complete model. We will moreover consider how to rewrite the gene associations during model inference.

Different kinds of models, with different degrees of sophistication, can be defined for the same organism. For this work we have selected four complementary use cases for model inference using gene associations. We will show that each requires a different kind of assertion on the reaction, and consequently a different kind of gene association. Consequently, the step (1 ) of figure 1, which projects a complete model to a minimal model annotated by gene associations, depends on the type of models we wish to infer.

• Model inference uses an annotated reference model and homology between the genes in that reference and the genes in a target organism, to infer a model for the target organism [cite]. The assertion that must be tested for each reaction is whether the genome of the target contains the required genes, for example those coding for the enzymes that catalyze the reaction.

• Flux balance analysis (FBA) calculates the fluxes through the reactions that maximize an objective function under a steady-state assumption [cite]. FBA can be used to simulate the effect of gene knockouts on metabolic fluxes. The assertion that must be tested is whether a reaction should be removed from the system because its requisite genes are knocked out or deactivated.

• Hierarchical modeling decomposes the complex cellular system into a hierarchy of independently-modeled stochastic transition systems, that communicate through flows and can synchronously perform state transitions through firing of global events [cite]. The assertion that must be tested is whether the guard on a transition, defined in terms of state variables such as protein expression or metabolite concentration, permits it to be executed.

• Dynamic simulation converts the metabolic model into a system of differential equations that define the evolution of molecule concentrations over time. Gene associations in dynamic simulation assert dependencies between variables that ultimately might be any molecular species involved in the model.

For each of these uses we define the appropriate gene associations and show how to derive them from a complete metabolic model. This derivation will provide a different gene association correctness criterion for each use case.

To illustrate the differences between the four use cases and the corresponding gene associations, we use three examples from Saccharomyces cerevisiae:

• Ribonuclease reductase, reaction RNR;

• Ferrocytochrome-C oxidoreductase, reaction CYC ;

• 3-methyl-2-oxopentanoate decarboxylase, reaction PDC.

Each of these examples presents a different challenge: the RNR system must take the cell cycle into account, the CYC system must take compartments into account, and the PDR system must take gene regulation into account.

By showing how the appropriate gene associations can be derived from a complete model, we provide the means to define a gold standard for annotation of a metabolic model. However, in practical applications we rarely have a complete model for which gene associations could be derived; rather, gene association annotations are defined at the same time as the definition or the inference of the biochemical reactions in the model. The interesting question in these applications is how to improve the quality of the reactions and their annotations. We show how the use-specific definition of gene associations provides a consistency check on the model being defined.

Definitions and Notation

Gene associations or Gene-Protein-Reaction (GPR) relations are boolean formulas constructed using the operators "AND" and "OR" over two levels: genes to proteins and proteins to reactions. In the first level, the formulas express how proteins are coded by their genes. In the second level, the formulas express the dependencies of reactions on proteins (Figure 2). Identifying GPR relations is essential for the reconstruction of metabolic networks[FP08, FST05, FHT + 09] and a critical step in that reconstruction.

Inria

The gene associations in a metabolic network can be identified through manual curation, a task that requires extensive experimental data [FGAT09] (e.g UniProtDB [C + 11]), or by automatic methods. The latter are founded on knowledge engineering, and are principally based on sequence similarity[KDWV13, WLT05, AJW + 08]. Because extensive experimental data are rarely available, efficient automatic methods are necessary to decipher correct annotations [RCO + 13].

Gene associations are stored in many public data banks such as KEGG [START_REF] Kanehisa | Kegg: kyoto encyclopedia of genes and genomes[END_REF] and MetaCyc [CFF + 06].

Association de gènes pour R : (G1 and G2) or (G3 and G4) The knowledge represented in the model by gene associations is indispensible for predicting cell phenotypes under genetic modifications such as gene knockouts in flux balance analysis (FBA) [START_REF] Paul A Jensen | Tiger: Toolbox for integrating genome-scale metabolic models, expression data, and transcriptional regulatory networks[END_REF], and for taking into account the effects of gene regulation in dynamic simulation [START_REF] Kumar | Metrxn: a knowledgebase of metabolites and reactions spanning metabolic models and databases[END_REF].

G4 G1 G2 G3 R OR G1 G2 G4 G3 & &

Gènes

Let us commence by defining the logical representation we will use to document a state of knowledge about a set of reactions in a metabolic model and a set of gene-protein-reaction relations associated with those reactions. In this representation, genes and proteins are atomic formulas (or atoms), and the relations between these objects and the state of the system are clauses. For example, the clause "A:B A, B" means that, to the complex A:B can only be formed if the proteins A and B are already present. We define inference rules that combine and rewrite clauses in order to derive the formulas that are required for each use.

The following notation is used for the assertions employed by the logical representation:

• (A, c) -the molecule A is present in compartment c
• A:B -a macromolecular complex is formed from A and B

• degrad(A) -the molecule A is degraded by the cell

• dis(A:B) -the complex A:B is dissociated by the cell

• ox(A) -the molecule A is oxidized • red(A) -the molecule A is reduced • ox-red(A) -the redox cycle of the molecule A is present
This representation permits one to write clauses of the kind C, D A, B, where A and B imply C and D, for example as a result of a reaction A + B → C + D. We will see, in what follows, that the formula required for each use can be obtained by rewriting clauses, up to a use-specific halting criterion.

BioRica [AGS + 11] is a high-level modeling framework that integrates discrete, continuous, stochastic, nondeterministic, and timed behaviors. It provides a means to hierarchically compose models in a nonambiguous way. BioRica extends the ideas of AltaRica [START_REF] Arnold | Construction and analysis of transition systems with MEC[END_REF] [APGR99] to biological systems, allowing specifications that use of differential equations as well as automata with stochastic and nondeterministic behaviors. BioRica is designed for modeling complex systems through the combination of existing models of subsystems, in particular those imported from SBML models.

The basic building block of BioRica is the node, which has discrete or continuous behavior. The specification of a node has eight fields: state, flow, sub, event, trans, assert, init and extern. The field state declares state variables and their value domains. Hierarchical connections to other BioRica nodea are specified using the field flow, which defines input-output connections, and the field sub, which lists nodes that the current node uses hierarchically. The field event declares the name of the discrete events that are possible, and the field trans defines the state transitions that are caused by these events. The field assert imposes constraints on input-output flows and on states. The field init declares the initial values of state variables. Finally, the field extern is used to specify event delay distributions and priorities between events.

Strategies for Deriving Gene Associations

In what follows we propose strategies for deriving gene associations for each of the use cases presented in 1. All of the strategies are based on a common method: substitution and simplification rules are applied top-down until a halting criterion is reached. The choice of halting criterion to use depends on the nature of the use case being studied.

In the case of model inference, to determine whether a reaction can occur in an organism, it suffices to know whether the enzymes that catalyze it are present or not. Consequently, gene association rewriting must halt as soon as the genes that code for the enzymes catalyzing the reaction are observed.

Steady-state flux balance analysis supposes that a reaction can carry any flux that is physiologically possible, as long as the necessary genes are expressed. Assuming steady-state, gene association reqriting must halt before any substitution, such as gene expression, for example, that would correspond to a change of state.

In the case of hierarchical modeling, the nodes in the BioRica specification correspond to a model that defines a natural limit of knowledge. Consequently, gene association rewriting is limited to the confines of each node. These limited can be identified syntactically in the model's specification.

Inria

Finally, in dynamic simulation, since any reaction may contribute to the simulation, gene association rewriting must continue as long as clauses can be rewritten, and can in principle concern the entire model. In the following, we will not illustrate this strategy.

The three examples that follow have been chosen to illustrate the process of deriving gene associations in metabolic networks. Each example starts from one central reaction of interest, but may include up to fifteen associated reactions, required for understanding the prerequisite conditions of the central reaction.

• RNR, the ribonucleotide reductase system modifying RNR by Trx1p, specifically the reaction catalyzed by the bisulfite Trx1p (r 8 dans Yeastnet).

• CYC, the complex Cyc7p:Ccp1p that catalyzes the transformation of ferricytochrome-C to ferrocytochrome-C (r 0469 dans Yeastnet).

• PDC, the transcription factor complex Thi3p:Pdc2p or Thi3p:Thi2p that induces the four paralogs of Pcd1p, who catalyze the transformation of 3-methyl-2-oxopentanoate en 2-methyl butanol ( r 0067 dans Yeastnet) 4 Ribonucleotide reductase system (RNR)

Ribonucleotide reductase (RNR) converts ribonucleotides to deoxiribonucleotides during the S phase of the cell cycle. Recall that the cell cycle corresponds to the set of steps constituting and delimiting the life of a cell, and is composed of several growth stages: the Mitose/Méiose phase, where the cell divides, the Gap1 phase, where the cell grows until it reaches a critical size that signals transition to the S phase, where DNA is replicated, before transitioning to the Gap2 phase.

During the Gap1, Gap2 and Mitose/Méiose phases, when ribonucleotide reduction is not necessary, the heterodimer complex Rnr2p:Rnr4p is imported into the nucleus through the action first of the addressing protein Dif1p, which brings the complex Rnr2p :Rnr4p in proximity to the nucleus, and second of the protein Kap122p that imports the complex into the nucleus. The Rnr2p:Rnr4p complex is retained in the nucleus by the anchoring protein Wtm1p.

During the S phase of the cell cycle, the Rnr2p:Rnr4p complex is moved from the nucleus to the cytoplasm by two cooperating mechanisms:

• The degradation Dif1p, which prevents entry of new Rnr2p:Rnr4p complexes into the nucleus.

• The action of the export protein Crm1p, which transports the complex Rnr2p:Rnr4p from the nucleus to the cytoplasm.

The Rnr2p:Rnr4p complex associates with the homodimère complex Rnr1p:Rnr1p to form the new complex Rnr1p:Rnr1p:Rnr2p:Rnr4p, which is then reduced by one of the other of the isozymes od thioredoxine Trx1p or Trx2p. These proteins are coded respectively by the TRX1 and TRX2. The reduced complex Rnr1p:Rnr1p:Rnr2p:Rnr4p in turn catalyzes the reduction of diphosphate ribonucleoside (NDP) into diphosphate deoxyribonucleoside (dNDP) and is thus oxidized.

A graphic representation of the model used here is presented in figure 3. Starting from this representation, we determined a logicial representation and a BioRica specification corresponding to the reaction NDP

∧ O 2 ∧ 2H 2 → dNDP ∧ 2H 2 O.
ox-red(Trx1p) ∨ ox-red(Trx2p) (Trx1p ∨ Trx2p) ∧ Trr1p (3)
Rnr1p:Rnr1p:Rnr2p:Rnr4p Rnr1p:Rnr1p ∧ Rnr2p:Rnr4p (5)

Rnr1p:Rnr1p Rnr1p ∧ Rnr1p (6)
Rnr2p:Rnr4p Rnr2p ∧ Rnr4p (7)

(Rnr1p:Rnr1p:Rnr2p:Rnr4p, c) (Rnr1p:Rnr1p, c) ∧ (Rnr2p:Rnr4p, c) (8)
(Rnr1p:Rnr1p, c) (Rnr1p:Rnr1p, n) ∧ Crm1p ∧ Phase(S) (9)
(Rnr1p:Rnr1p, n) (Rnr1p:Rnr1p, c) ∧ Kap122p ∧ Dif1p ∧ ¬Phase(S) (10)

BioRica Description

The hierarchical model of the RNR system is comprised of three BioRica nodes:

• Subcomplex, which describes the formation of the Rnr2p:Rnr4p heterodimer and its importation into the nucleus. We suppose that the formation of the complex is the limiting step in the reaction; the kinetics of the addressing and importing reactions are assumed to be much more rapid, and are not explicitly modeled.

• Complex, which describes the formation of the Rnr1p:Rnr1p:Rnr2p:Rnr4p heterotetramer in the cytoplasm during S phase. We suppose that this formation is only limited by the presence of Rnr2p:Rnr4p in the cytoplasm, which is in turn controlled by the transition to S phase. The transporter Crm1 is not explicitly modeled.

• Reaction, which described the system of reactions neighboring ribonucleotide reductase. This includes the cycles of reduction-oxidation of the reductase heterotetramer and both of the thioredoxines Trx1 and Trx2.

The Main node links these three sub-nodes by boolean variables that indicate the presence in the cytoplasm of the heterodimer Rnr2p:Rnr4p and the heterotetramer Rnr1p:Rnr1p:Rnr2p:Rnr4p.

It also synchronizes the transitions to the S phase by the sub-nodes Subcomplex and Complex.

The complete BioRica specification is as follows.

) + ln( C(O 2 ) ic(O 2 ) ) -2 ln( C(H 2 O) ic(H 2 O) ) -ln( C(dN DP ) ic(dN DP ) ))); trans Inria

Gene Associations for Model Inference

We can conclude that the ribonucleotide reductase reaction can occur in a target organism if there is sufficient evidence to rewrite its gene association. The reaction is catalyzed by the oxidation cycle of the Rnr1p:Rnr1p:Rnr2p:Rnr4p complex, so in order in instantiate the reaction we must be assured that the genes coding for the constituant proteins of the complex are present in the target organism (figure 4). The gene association of the ribonucleotide reductase reaction is given by RNR1 ∧ RNR2 ∧ RNR4. This is derived formally from the presence of the oxidation cycle of the Rnr1p:Rnr1p:Rnr2p:Rnr4p complex: ox-red(Rnr1p:Rnr1p:Rnr2p:Rnr4p)

Rnr1p:Rnr1p:Rnr2p:Rnr4p ∧ (ox-red(Trx1p) ∨ ox-red(Trx2p)) eq(3) Rnr1p:Rnr1p ∧ Rnr2p:Rnr4p eq(5) Rnr1p ∧ Rnr2p ∧ Rnr4p eqs(6) (7) RNR1 ∧ RNR2 ∧ RNR4 eqs(18) (19) (20)

Gene associations for Flux Balance Analysis

Flux Balance Analysis takes into account fluxes through reactions, under the hypothesis that no state change occurs. In consequence, we may not apply substitution or simplification rules that lead to gene expression (e.g. eq(22)) or to formation or dissociation of protein complexes (e.g. eq (5)). In this precise case, we suppose that the enzymes and enzymatic complexes are present, and we do not attempt to take into account the cycles of reduction reactions (Figure 5). The gene association of the ribonucleotide reductase reaction is given by (Trx1p ∨ Trx2p) ∧ Trr1p. This is derived formally from the presence of the oxidation cycle of the Rnr1p:Rnr1p:Rnr2p:Rnr4p complex using rules that do not compise gene expression or protein complex formation: ox-red(Rnr1p:Rnr1p:Rnr2p:Rnr4p) Rnr1p:Rnr1p:Rnr2p:Rnr4p ∧ ox-red(Trx1p) ∨ ox-red(Trx2p) eq(3) Rnr1p:Rnr1p:Rnr2p:Rnr4p ∧ (Trx1p ∨ Trx2p) ∧ Trr1p eq(4)

Gene Associations for Hierarchical Modeling

The reduction of ribonucleotide to deoxyribonucleotide is performed by the Rnr1p:Rnr1p:Rnr2p:Rnr4p heretotetramer complex, which is formed by the association of the Rnr1p:Rnr1p homodimer with the Rnr2p:Rnr4p heterodimer released from the nucleus. The reductase complex is activated by a Trx1p or Trx2p thioredoxine, which is in turn activated by the Trr1p thioredoxine reductase. These reactions occur at the same time: each ribonucleotide reduction exidizes the reductase complex, which is reduced by the thioredoxine, itself reduced by the Inria thioredoxine reductase. Thus, to describe these dependencies, we require the genes coding the Rnr1p:Rnr1p:Rnr2p:Rnr4p heterotetramer reductase complex, the thioredoxine that performs the reduction, as well as the thioredoxine reductase (Figure 6). The Cyc1p and Cyc7p isoforms of cytochrome c are electron transport proteins taking part in the electron transport chain in mitochondria. In this system, the Cyc1p or Cyc7p proteins form a complex with cytochrome c peroxidase Ccp1p. In this complex, Ccp1p c catalyzes the oxidation of 2-ferrocytochrome c, which is the reduced form of Cyc1p or Cyc7p (figure 7). Les cytochrome c isoformes Cyc1p et Cyc7p sont des protéines de transport d'électrons dans la chaîne de transport d'électrons mitochondriale. Dans ce système, la protéine Cyc1p ou Cyc7p forme d'abord un complexe avec cytochrome c peroxidase Ccp1p, et dans ce complexe, Ccp1p catalyse l'oxydation du 2-ferrocytochrome c, qui est la forme réduite de Cyc1p or Cyc7p (Figure 7).

Gene Associations for Model Inference

The genes needed so that this reaction may occur in the target organism are the genes coding for the catalyst Cc1p and the complex Cyc1p|Cyc7p:Ccp1p (figure 8). Ccp1p ∧ Cyc7p:Ccp1p eq(29) Ccp1p ∧ Cyc7p eq(30) CCP1 ∧ CYC7 eqs(32) (31)

Gene Associations for Flux Balance Analysis

In this example no state changes occur, so the reaction may be instantiated if the enzyme and the complex are already present. (figure 9). R Ccp1p ∧ Cyc1p:Ccp1p eq(25) R'

Ccp1p ∧ Cyc7p:Ccp1p eq(29)

Gene Associations for Hierarchical Modeling

The oxidation of cytochrome c (Cyc1p ou Cyc7p) by Ccp1p takes place within the Cyc1p:Ccp1p or the Cyc7p:Ccp1p) complex, which associates cytochrome c with cytochrome c peroxidase (figure 10). R Ccp1p ∧ Cyc1p:Ccp1p eq(25) Ccp1p ∧ Cyc1p eq(26) CCP1 ∧ CYC1 eqs(28) (27) The enzymes of 3-methyl-2-oxopentanoate decarboxylase, Pdc1p, Pdc5p, Pdc6p, et ARO10, react to 2-oxo-acids like 3-methyl-2-oxopentanoate, which is an intermediate in the synthesis of isoleucine, (figure 11).

Thi3p is a regulatory protein of this reaction. It binds to the transcription factors Thi2p and Pdc2p. In the presence of thiamine, the Pdc2p:Thi3p and Thi2P:Thi3p complexes dissociate, and thus deactivate the transcription factors. The genes PDC1 and PDC5, which code respectively for the proteins Pdc1p and Pdc5p, are regulated by the Pdc2p:Thia3p complex or by the Thi2p:Thi3p complex, and by the absence of the protein Pdc1p. However, the genes PDC6 and ARO10, which code respectively for the proteins Pdc6p and Aro10p, are not regulated by these mechanisms. Pdc1p ∨ Pdc5p ∨ Pdc6p ∨ Aro10p eq(33)

Gene Associations for Hierarchical Modeling

The transformation of 3-méthyl-2-oxopentanoate into 2-méthyl butanol is performed by one of the enzymes Pdc1p, Pdc5p, Pdc6p, or Aro10p (figure 14) .

R(PDC)

Pdc1p ∨ Pdc5p ∨ Pdc6p ∨ Aro10p eq(33) PDC1 ∨ PDC5 ∨ PDC6 ∨ ARO10 eqs(34),(35),(40),(41)

Discussion

This report defines a formal criterion for correctness of gene associations: a gene association is correct if it is the one that would have been derived from a complete model. Our goal has been to define this criterion, and not so much to provide a method for deriving the gene association. In most practical cases a complete model is not available.

Starting from the definition of gene associations, the requirements of specific use cases, and the possible influences of a (gene or protein) variable, we defined simplification and substitution rules that produce the most pertinent gene association for a given reaction. These rules may also be used to provide guidelines for curators who build gene associations by hand.

The criterion may be used to define methods for finding gene associations through the analysis of experimental results. It may also be used to evaluate the form of boolean formulas used in gene associations, in order to facilitate knowledge extraction by inference procedures.

Figure 2 :

 2 Figure2: (left) The gene association for a reaction R catalyzed by a protein complex comprised either of the aggregation of products of genes G1 and G2, or of the aggregation of the products of G3 and G4. (right) The formula for the gene association, written as a binary tree.

Figure 3 :

 3 Figure 3: Explicit model for the ribonucleotide reductase reaction (RNR)

  ) → ox(Trx1p | Trx2p) := true; (Rnr1p:Rnr1p:Rnr2p:Rnr4p, c),ox(Trx1p | Trx2p) ox(RNR) → ox(Rnr1p:Rnr1p:Rnr2p:Rnr4p, c):= true; ox(Rnr1p:Rnr1p:Rnr2p:Rnr4p, assert C.(Rnr1p:Rnr1p:Rnr2p:Rnr4p, c) = R.(Rnr1p:Rnr1p:Rnr2p:Rnr4p, c); N.(Rnr2p:Rnr4p, c) = C.(Rnr2p:Rnr4p, c); sync <N.phaseS, C.phaseS>; edon

Figure 4 :

 4 Figure 4: Explicit model of the ribonucleotide reductase reaction, outlining in red the part of the model considered when performing model inference

Figure 5 :

 5 Figure 5: Explicit model of the ribonucleotide reductase reaction, outlining in red the part of the model considered when performing flux balance analysis

Figure 6 :

 6 Figure 6: Explicit model of the ribonucleotide reductase reaction, outlining in red the part of the model considered when performing hierarchical modeling

Figure 7 :

 7 Figure 7: Explicit model of the ferrocytochrome c reaction

Figure 8 :

 8 Figure 8: Explicit model for the Ferrocytochrome c system, outlining in red the part of the model considered when performing model inference R Ccp1p ∧ Cyc1p:Ccp1p eq(25) Ccp1p ∧ Cyc1p eq(26) CCP1 ∧ CYC1 eqs(28) (27) R' Ccp1p ∧ Cyc7p:Ccp1p eq(29) Ccp1p ∧ Cyc7p eq(30) CCP1 ∧ CYC7 eqs(32) (31)

Figure 9 :

 9 Figure 9: Explicit model of the Ferrocytochrome c system, outlining in red the part of the model considered when performing flux balance analysis

Figure 11 :Figure 12 :

 1112 Figure 11: Explicit model of the 3-méthyl-2-oxopentanoate décarboxylase reaction

Figure 13 :

 13 Figure 13: Explicit model of the 3-methyl-2-oxopentanoate decarboxylase system, outlining in red the part of the model considered when performing flux balance analysis

RR n°9246

Logical Description

Ccp1p CCP1 (28)

Ccp1p CCP1 (32)

BioRica Description

The hierarchical model of the CYC system contains a single hybrid node Main, composed of discrete transitions that represent gene expression and the formation of the Ccp1p:cytochrome c complex, and one continuous equation describing the oxidation of cytochrome c by Ccp1. The oxidation of Ccp1 may only occur once the Ccp1p:cytochrome c complex is formed.

BioRica Description

The hierarchical model of the PDC system is composed of two sub=nodes:

• Complex, which describes the association of Thi3p with either Thi2p or with Pdc2p, to form the transcription factor that regulates the expression of the genes PDC1 and PDC5.

• Factor, which is a hybrid node composed of discrete transitions that represent the expression of the genes PDC1, PDC5, PDC6, or ARO10, and a continuous equation representing the decarboxylase reaction. The expression of genes PDC1 and PDC5 depends on the transcription factor.

Gene Associations for Model Inference

The reaction is catalyzed by one of the enzymes Pdc1p, Pdc5p, Pdc6p, and Aro10p, which are products of the genes PDC1, PDC5, PDC6, and ARO10. Consequently, the reaction may be inferred if at least one of these genes is present (figure 12). R(PDC) Pdc1p ∨ Pdc5p ∨ Pdc6p ∨ Aro10p eq(33) PDC1 ∨ PDC5 ∨ PDC6 ∨ ARO10 eqs(34) (35) (40) (41)

Gene Associations for Flux Balance Analysis

In this specific case, we suppose that the enzymes and enzymatic complexes are present (or not), and we only need to model the cycles of reduction reactions.

Inria