Riccardo Porreca
email: riccardo.porreca@control.ee.ethz.ch

Eugenio Cinquemani

John Lygeros

Giancarlo Ferrari-Trecate

Invalidation of the structure of genetic network dynamics: A geometric approach

Keywords: systems biology, identification, quasiconvexity, unate functions, sigmoidal activation functions

This work concerns the identification of the structure of a genetic network model from measurements of gene product concentrations and synthesis rates. In earlier work, we developed a data preprocessing algorithm that is able to reject many hypotheses on the network structure by testing certain monotonicity properties for a wide family of network models. Here we develop a geometric interpretation of the method. Then, for a relevant subclass of genetic network models, we extend our approach to the combined testing of monotonicity and convexity-like properties associated with the network structures. The theoretical aspects and practical performance of the enhanced methods are illustrated by way of numerical results.

INTRODUCTION

Genetic networks govern the behavior of living cells in response to changes in the environment, and determine growth, replication, and death of cells. They are composed of genes, i.e. pieces of the DNA strand encoding a specific protein. Proteins are synthesized in several copies upon gene expression and participate in the regulation of the expression of other genes, thus giving rise to a complex network of biochemical interactions.

Modern technologies for the time-course measurement of gene expression, such as gene reporter systems, allow one to step from pure topological modelling of gene networks to the modelling of the interaction dynamics. However, this requires setting up a dynamical model of the network whose structure and parameters are typically unknown or uncertain. Data-based identification of an accurate model is challenging due to the size of the family of possible model alternatives. Yet, a priori biological knowledge may be exploited so as to ameliorate the complexity of the problem.

In [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF], we developed an identification strategy for genetic network dynamics with a unate structure. These are ODE models built upon a family of Boolean network models which reportedly capture most of the experimentally observable gene regulation interactions [START_REF] Grefenstette | An analysis of the class of gene regulatory functions implied by a biochemical model[END_REF]. In [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF], we showed that unate models possess monotonicity properties that can be tested inexpensively on experimental data, so as to discard entire sets of model hypotheses and focus the search on model structures consistent with the data. A similar approach to the identification of gene network topology using time-course data was recently developed in [START_REF] Julius | Genetic regulatory network identification using monotone functions decomposition[END_REF], relying on discrete-time models still having monotonicity properties and without taking explicitly into account noise affecting the data. Other relevant approaches to reverse engineering of gene networks, based on different modeling assumptions and inference techniques, can be found in [START_REF] Bansal | How to infer gene networks from expression profiles[END_REF][START_REF] Cantone | A yeast synthetic network for in vivo assessment of reverse-engineering and modeling approaches[END_REF][START_REF] Anderson | On validation and invalidation of biological models[END_REF][START_REF] Julius | Genetic network identification using convex programming[END_REF][START_REF] August | Efficient, sparse biological network determination[END_REF][START_REF] Zavlanos | Inferring stable genetic networks from steady-state data[END_REF].

One question that arises naturally is whether additional properties of the models in the class (other than monotonicity) can be exploited so as to further narrow down the search for valid models. In this paper we address this question by considering a subclass of unate models. In particular, we leverage on the analysis in [START_REF] Nikolajewa | Boolean networks with biologically relevant rules show ordered behavior[END_REF] showing that 87% of the 139 Boolean genes regulatory rules considered in [START_REF] Harris | A model of transcriptional regulatory networks based on biases in the observed regulation rules[END_REF] belong to a narrower class, which we refer to as S 0 ∪ S 1 , and show that ODE models with S 0 ∪ S 1 structure possess convexity-like properties that can be used for checking the consistency of different model hypotheses with the experimental data. To this purpose, we introduce a geometric framework that also provides an alternative interpretation of the methods by [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF]. In the economics literature, a similar approach was considered by [START_REF] Hanoch | Testing the assumptions of production theory: A nonparametric approach[END_REF] for testing hypotheses on production processes.

In Sections 2.1 and 2.2 we review our results from [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] on modelling and invalidation of genetic network dynamics with unate structure. In Section 2.3, we give these results a new geometrical interpretation. The same approach is exploited in Section 3 to analyze convexity-like properties of the models with S 0 ∪ S 1 structure and set up new model invalidation strategies. In Section 4 we discuss efficient implementations of these methods as well as a strategy to handle noisy data. Theoretical and experimental results are discussed in Section 5 by way of illustrative simulations. Mathematical proofs of the results developed in the paper can be found in [START_REF] Porreca | Invalidation of the structure of genetic network dynamics: A geometric approach[END_REF]. Sections 2.3 and 3-5 provide the original contributions of this work.

Notation: If v ∈ R d , diag(v) stands for the diagonal matrix V ∈ R d×d with V ii = v i , i = 1, . . . , d. For a set M ⊂ R n , Conv(M ) denotes its convex hull. The cardinality of a finite set P is |P |.

UNATE MODELS AND MONOTONICITY

Sections 2.1 and 2.2 are a concise review of our results in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] and provide the context and the background for the results developed in this work. Section 2.3 reconsiders these results from a geometrical point of view and establishes the ground for the results presented in the later sections.

Genetic network models with unate structure

In the context of Boolean network modelling, the activation status of gene i, with i = 1, . . . , n, in a network of n genes is encoded by a binary variable X i that takes the value 1 if the gene is active and 0 otherwise. A Boolean rule B i (X) : {0, 1} n → {0, 1} describes the logics governing the activation of gene i as a function of X = (X 1 , . . . , X n ). Depending on what genes regulate the expression of gene i, B i will effectively depend on only some of the entries of X. Based on biochemical reaction modelling arguments [START_REF] Grefenstette | An analysis of the class of gene regulatory functions implied by a biochemical model[END_REF], it can be argued that most regulatory interactions are well described by unate functions, i.e. Boolean functions that are monotone (either nondecreasing or nonincreasing) in each of the input variables, meaning that increasing the value of a single X i from 0 to 1 can only force a specific trend of B i (X) (either increasing from 0 to 1 or decreasing from 1 to 0) or leave it unchanged, regardless of the value of the other variables. In conjunctive normal form, unate functions are given by

B i (X) = hi l=1 j∈J il Xj , (1)
where "∧" and "∨" stand for "and" and "or", respectively, h i is a nonnegative integer and each J il is a nonempty set of indices from {1, . . . , n}. Furthermore, unate functions must fulfill the additional constraint that each variable Xj appears in (1) exclusively as either X j or ¬X j , where "¬" is negation. By convention, a conjunction of h i = 0 terms is equal to 1. An example of function that is not in this class is the exclusive or ("xor"). Indeed, if X 2 = 0 then X 1 xor X 2 increases from 0 to 1 as X 1 goes from 0 to 1, while if X 2 = 1 the function decreases from 1 to 0. Moreover, the conjunctive normal form

X 1 xor X 2 = (X 1 ∨ X 2 ) ∧ (¬X 1 ∨ ¬X 2
) includes e.g. both X 1 and its negated form ¬X 1 . The class of unate models represents all networks where each regulator of a given gene acts unambiguously either as an activator or as a repressor of that gene, though it may promote expression of one gene and inhibit the expression of another gene. Unate functions include the so-called Hierarchically Canalizing Functions (HCFs) [START_REF] Aracena | Maximum number of fixed points in regulatory boolean networks[END_REF], which capture a large class of the known regulatory interactions among genes and are intimately related with the stability properties of the network [START_REF] Kauffman | Genetic networks with canalyzing boolean rules are always stable[END_REF][START_REF] Szallasi | Modeling the normal and neoplastic cell cycle with "realistic boolean genetic networks": their application for understanding carcinogenesis and assessing therapeutic strategies[END_REF].

In several cases of interest, further assumptions can be made a priori on the structure of B i (X).

Example 1

It is shown in [START_REF] Nikolajewa | Boolean networks with biologically relevant rules show ordered behavior[END_REF] that many gene activation rules fall into two subclasses of HCF, S 0 and S 1 , having the following form:

B i (X) = Xj1 ∧ Xj2 ∧ Xj3 ∧ • • • ∧ Xj ℓ , if B i ∈ S 0 , [ Xj1 ∨ Xj2 ] ∧ Xj3 ∧ • • • ∧ Xj ℓ , if B i ∈ S 1 , ( 2
)
where ℓ is the number of effective inputs of B i (X) and j 1 , . . . , j ℓ are pairwise different indices from the set {1, . . . , n}.

We are interested in quantitative models of gene expression. We consider models of the form [START_REF] De Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF]

ẋi = g i (x) -γ i (x) , (3)
where i = 1, . . . , n denotes the ith of n genes, x i ≥ 0 denotes the concentration of the corresponding product, x = (x 1 , . . . , x n ), and g i (x) ≥ 0 and γ i (x) ≥ 0 are synthesis and degradation rates, respectively. Functions g i (x) and γ i (x) are generally used to model the regulatory effects on the synthesis and degradation of the ith gene of the network. In view of the falsification approach developed in this work, we are especially interested in synthesis rate regulation functions. In particular, for i = 1, . . . , n, we assume that

g i (x) = κ 0,i + κ 1,i b i (x) , (4)
where κ 0,i and κ 1,i are nonnegative constants and the gene activation level b i (x) is typically a combination of switch-like (e.g. sigmoidal) regulatory functions describing the effect of protein j on the expression of gene i and the synthesis of the corresponding protein. In order to account for the discrete regulatory logics (1) in the quantitative model ( 3)-( 4), we follow an approach inspired by [START_REF] Plahte | A methodological basis for description and analysis of systems with complex switchlike interactions[END_REF]. Each variable X i is replaced by a monotone, nondecreasing sigmoidal function σ + : [0, +∞) → [0, 1] of the concentration x i . Given any two functions τ (x) and τ ′ (x) representing the Boolean expressions T (X) and T ′ (X), ¬T (X) is replaced by 1τ (x) and T (X) ∧ T ′ (X) by τ (x) • τ ′ (x). In particular ¬X i is represented by σ -(x i ) = 1σ + (x i ). Applying these rules to (1) leads to the following class of models.

Definition 1

A unate model is given by ( 3) and ( 4) where, for some integer h i ≥ 0 and some sets of indices

J il ⊆ {1, . . . , n}, with l = 1, . . . , h i , b i (x) = hi l=1 1 - j∈J il 1 -σ ± (x j ) , (5)
where either σ ± (x j ) = σ + (x j ) or σ ± (x j ) = σ -(x j ) and, by convention, products over an empty set return 1.

A typical choice of sigmoid is the Hill function [START_REF] Yang | An analytical rate expression for the kinetics of gene transcription mediated by dimeric transcription factors[END_REF][START_REF] Keller | Model genetic circuits encoding autoregulatory transcription factors[END_REF]. For any choice of the cooperativity parameter d ≥ 1 and the threshold parameter η ≥ 0, this function increases monotonically from 0 to 1, satisfies σ + (η) = 1/2 and dσ + dx (η) ≥ 0 increases with d. For d = 1, in particular, one recovers Michaelis-Menten kinetics. An alternative choice is the logistic function

σ + (x) = 1/[1 + (η/x) d ]
σ + (x ′ ) = 1 1 + e -d(x ′ -η) , x ′ ∈ R ≥0 , (6)
with analogous interpretation for d and η. Note that in this case σ + (x ′ ) = 0 when x ′ = 0. If one looks at (6) over the entire R, it is immediate to see that Hill functions are in one-to-one correspondence with (6) via the transformations

x ′ = log(x) , (7)
η = log(η) . (8)
Therefore, we accept that b i and g i are generally defined over the entire R n and assume without loss of generality that all sigmoids in (5) are logistic functions. The methods developed below also apply to models involving Hill functions, by taking logarithms of concentration variables and measurements, as above.

Invalidation of unate models: sign patterns

Consider the problem of identifying a unate model for gene i from the dataset D = {(x k , g k i ) : k = 1, . . . , m}, where each x k is a vector of protein concentrations and g k i = g i (x k ) is the corresponding synthesis rate. In practice, (noisy versions of) measurements D can be obtained by perturbation experiments (see [START_REF] Sontag | Network reconstruction based on steady-state data[END_REF] and references therein) or time-course experiments [START_REF] De Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]. In particular, methods proposed in this paper are well adapted to gene reporter systems where average promoter activities over a cell population are sampled with relatively high frequency. In fact, protein concentrations x and synthesis rates g i can be inferred from coarse promoter activity data e.g. by means of the nonparametric estimation methods proposed in [START_REF] De Jong | Experimental and computational validation of models of fluorescent and luminescent reporter genes in bacteria[END_REF]. Since x and g i are both observed, one faces a regression problem for the function g i only, i.e. the specific form of γ i (x) in (3) is irrelevant.

A fundamental source of complexity is that the function b i (x) in (5) depends upon discrete quantities, i.e. the integers h i , the sets J il and the signs of the sigmoids. For realistic size problems, it is computationally prohibitive to search for the best fitting model by identifying values for all the parameters (κ 0,i , κ 1,i , thresholds, cooperativity parameters) for all possible combinations of discrete quantities. Hence, we focus on the problem of invalidating families of unate models on the basis of the dataset D independently of the value of continuous parameters.

For ease of exposition, we start by assuming that data are noiseless and return to the case where measurement noise is present in Section 4.1. Moreover, since the problem is the same for all genes, we drop the index i to simplify the notation. In [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] we addressed the invalidation problem by exploiting monotonicity properties of g(x). Given a model in the form (4)-( 5), we define its sign pattern p = (p 1 , . . . , p n ) ∈ {-1, 0, 1} n by

p j =      0, if j / ∈ J l , l = 1, . . . , h, 1, if σ ± (x j ) = σ + (x j ), -1, if σ ± (x j ) = σ -(x j ).
Note that several alternative structures of (5) share the same sign pattern p, e.g. b

(x) = σ + (x 1 )σ -(x 2 ) and b(x) = 1 -(1 -σ + (x 1 )) (1 -σ -(x 2 )
), corresponding to the Boolean formulas X 1 ∧ ¬X 2 and X 1 ∨ ¬X 2 , respectively, are both represented by p = (1, -1). For this reason, we introduce the family U (p) of unate models g(x) given by ( 4)-( 5) with sign pattern p. In light of Definition 1, g ∈ U (p) (and, similarly, the corresponding function b) is nondecreasing (respectively, nonincreasing) in x j if p j = 1 (respectively, p j = -1), and is independent of x j if p j = 0.

Proposition 1

For any g ∈ U (p) and l, k ∈ {1, . . . , m} it holds that

p j (x k j -x l j ) ≥ 0, j = 1, . . . , n ⇒ g(x k ) -g(x l ) ≥ 0
where x k j indicates the jth entry of x k . Based on Proposition 1, the family U (p) is invalidated by D if two indices l, k ∈ {1, . . . , m} exist such that p j (x k jx l j ) ≥ 0, j = 1, . . . , n and g kg l < 0 .

In this case, p is said to be inconsistent (with D).

Remark 1

In [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] we further introduced the concept of subpattern, i.e. p ′ is a subpattern of p (and p is a superpattern of p ′ ) if all nonzero entries of p ′ are equal to the corresponding entries of p. We showed that if p is inconsistent (respectively, consistent) then every subpattern (respectively, superpattern) p ′ is also inconsistent (respectively, consistent). This partial ordering relationship allowed us to characterize entire hierarchies of consistent patterns by way of minimal elements, and devise efficient algorithms for the storage and the enumeration of all consistent patterns.

A geometric approach to the invalidation of unate models

For a real-valued function g and ε ∈ R, define the super-level set T ε (g) = {x : g(x) ≥ ε} and the sub-level set B ε (g) = {x : g(x) ≤ ε}. We will now show that testing if a family U (p) is invalidated by D can be done using inner approximations of sets T g k (g), k ∈ {1, . . . , m}, computed on the basis of the dataset D only. For x ∈ R n and p ∈ {-1, 0, 1} n , define the cone

≥ (x, p) = {z ∈ R n : p j z j ≥ p j x j , j = 1, . . . , n}
with vertex x (see Fig. 1), that is the orthant defined by the nonzero entries of p, translated to x.

Intuitively, for any g ∈ U (p), p determines the direction of growth of g, hence ≥ (x k , p) is a region where g must be no smaller than g k . For a nonempty set D ′ ⊆ D, let K(D ′ ) ⊆ {1, . . . , m} be the set of indices of the elements of D contained in D ′ and define

M (D ′ , p) = k∈K(D ′ ) ≥ (x k , p) , µ(D ′ ) = min k∈K(D ′ ) g k .
The following result, that is illustrated in Fig. 1, shows that M (D ′ , p) provides a data-based inner approximation of the set T µ(D ′ ) (g).

Proposition 2

For all nonempty sets

D ′ ⊆ D, if g ∈ U (p) then M (D ′ , p) ⊆ T µ(D ′ ) (g).
It is important to notice that the set M (D ′ , p) depends on p but not on the particular g ∈ U (p). This allows us to redefine inconsistent sign patterns as follows.

Definition 2 A sign pattern p is m-inconsistent † if there is a nonempty set D ′ ⊆ D and (x * , g * ) ∈ D \ D ′ such that x * ∈ M (D ′ , p) and g * < µ(D ′ ). Otherwise p is m-consistent.
It can be shown that Definition 2 is equivalent to the definition of inconsistent sign pattern of Section 2, which means that one can only consider singleton sets D ′ . In particular, the construction of the falsification region M (D ′ , p), relies only on the monotonicity properties of g. However, the added benefit of this geometric approach becomes apparent when one focuses on subclasses of unate models, where larger falsification regions can be constructed using additional model properties. This is the approach we will follow in the next section.

QUASI-CONVEXITY ANALYSIS OF GENETIC NETWORK MODELS

Following on Example 1, we introduce subsets of unate models that are the algebraic counterpart of the S 0 ∪ S 1 Boolean models.

Definition 3

An S 0 model is given by ( 3) and ( 4) where, for some ℓ ∈ {0, . . . , n} and some subset {j

1 , . . . , j ℓ } of {1, . . . , n}, b(x) = σ ± (x j1 )σ ± (x j2 )σ ± (x j3 ) • • • σ ± (x j ℓ ) . (10)
Similarly, an

S 1 model is characterized by b(x) = b ∨ (x)b ∧ (x) , ( 11a
) b ∨ (x) = 1 -1 -σ ± (x j1 ) 1 -σ ± (x j2 ) , ( 11b
) b ∧ (x) = σ ± (x j3 )σ ± (x j4 ) . . . σ ± (x j ℓ ) . ( 11c
)
Finally an S 0 ∪ S 1 model is given by ( 3) and (4) if either [START_REF] Nikolajewa | Boolean networks with biologically relevant rules show ordered behavior[END_REF] or [START_REF] Harris | A model of transcriptional regulatory networks based on biases in the observed regulation rules[END_REF] holds.

In the sequel, S 0 (p) will denote the family of S 0 models g(x) given by ( 4) and [START_REF] Nikolajewa | Boolean networks with biologically relevant rules show ordered behavior[END_REF] with sign pattern p. Note that p defines an S 0 model up to the values of the kinetic and sigmoid parameters. In the case of S 1 models, the structure is parametrized by triplets (j ∨ , j ∧ , p), where p is a sign pattern and j ∨ = {j 1 , j 2 }, j ∧ = {j 3 , . . . , j ℓ } are sets of indices partitioning {i ∈ {1, . . . , n} : p i = 0}. We denote by S 1 (j ∨ , j ∧ , p) the family of S 1 models sharing the same structure (j ∨ , j ∧ , p), while S 1 (p) is the union of all families S 1 (j ∨ , j ∧ , p) sharing the same sign pattern p.

Remark 2

Note that if we allow for j ∨ = ∅ (and j ∧ = {i ∈ {1, . . . , n} : p i = 0}) in S 1 models, we recover S 0 models. This suggests that one can focus on the invalidation of families of S 1 models only. However, invalidation methods for S 0 models are simpler and hence, for the sake of clarity, we will discuss the S 0 and S 1 cases separately.

Next, we will show that S 0 ∪ S 1 models have quasi-convexity properties that can be used for invalidating entire families S 0 (p) or S 1 (j ∨ , j ∧ , p) using data in D only. Let us begin with some basic definitions and results of convex analysis (see [START_REF] Boyd | Convex Optimization[END_REF] for more details). i. for every α ∈ [0, 1] and every x, y ∈ D,

Definition 4

Let D ⊆ R n be a convex set. A function g : D → R is quasi-convex if
g αx + (1 -α)y ≤ max{g(x), g(y)} ; (12)
ii. for every ε ∈ R, the sub-level set B ε (g) is convex.

g is quasi-concave if -g is quasi-convex, that is: i'.
for every α ∈ [0, 1] and every x, y ∈ D,

f αx + (1 -α)y ≥ min{g(x), g(y)} ; (13)
ii'. for every ε ∈ R, the super-level set T ε (g) is convex.

The following proposition is the basis for new invalidation procedures.

Proposition 3

Function b(x) in ( 10) is quasi-concave (with respect to (x j1 . . . , x j ℓ )). Function b ∨ (x) in (11b) is quasi-convex (with respect to (x j1 , x j2 )), while b ∧ (x) in (11c) is quasi-concave (with respect to (x j3 . . . , x j ℓ )).

Note that Proposition 3 strongly relies on the standing assumption that sigmoids are logistic functions (see Section 2) but is completely independent of the signs of the sigmoids and the values of the cooperativity and threshold parameters. Moreover, quasi-convexity is not affected by multiplication by and addition of nonnegative scalars and hence g(x) is quasi-convex if and only if b(x) has the same property. In practice, this will allow us to infer properties of b(x) from data generated by the function g(x). We now apply these results to the invalidation of S 0 models (Section 3.1) and S 1 models (Section 3.2).

Invalidation of S 0 models

For any D ′ ⊆ D and any sign pattern p let C(D ′ , p) = Conv(M (D ′ , p)). The following result, that is illustrated in Fig. 2, shows that C(D ′ , p) provides a data-based inner approximation of the set T µ(D ′ ) (g) and hence leads to a new definition of inconsistent S 0 models.

Proposition 4 If g ∈ S 0 (p), then C(D ′ , p) ⊆ T µ(D ′ ) (g). Definition 5 A sign pattern p is c-inconsistent ‡ if there is a (nonempty) set D ′ ⊆ D and (x * , g * ) ∈ D \ D ′ such that x * ∈ C(D ′ , p) and g * < µ(D ′ ). Otherwise p is c-consistent.
Definition 5 strengthens Definition 2 for S 0 models since C(D ′ , p) ⊇ M (D ′ , p). In other words, if p is m-inconsistent, no g ∈ U (p) (and hence no g ∈ S 0 (p)) can generate the dataset D. However, if p is m-consistent the whole family S 0 (p) is still invalidated if p is c-inconsistent. Relations between m-and c-inconsistency are illustrated in the following example.

Example 2

With reference to Fig. 2, consider the dataset D = D ′ ∪ {(x # , g * )}, D ′ = {(x i , g i ), i = 1, . . . , 4} and the sign pattern p = (-1, -1). If g * < min i=1,...,4 g i one has x # ∈ M (D ′ , p) and g * < µ(D ′ ) and therefore, according to Definition 2, the sign pattern p is m-inconsistent. This means that all models in U (p) must be rejected, including

k 0 + k 1 σ -(x 1 )σ -(x 2 ) that is the only model in S 0 (p). If instead D = D ′ ∪ {(x * , g * )} because x * ∈ C(D ′ , p
), according to Definition 5, the sign pattern p is c-inconsistent. This means that the model in S 0 (p) has to be rejected. However x * ∈ M (D ′ , p) and therefore no model in U (p) can be invalidated based on m-inconsistency of p.

Invalidation of S 1 models

For S 1 models a convexity-like property does not globally hold. Hence the goal is to combine the different properties of (11b) and of (11c). There are different ways to do so, each leading to different conditions for the invalidation of model structures. For a generic z ∈ R d and p ∈ {-1, 0, 1} d , recall the definition of the cone

≥ (z, p) = {z ′ ∈ R d : p j z ′ j ≥ p j z j , ∀j = 1, . . . , d} (14a)
and define the cone

≤ (z, p) = ≥ (z, -p) = {z ′ ∈ R d : p j z ′ j ≤ p j z j , ∀j = 1, . . . , d} . (14b)
For sets of indices j ∨ = {j 1 , j 2 }, j ∧ = {j 3 , . . . , j ℓ } and a sign pattern p, let p ∨ = (p j1 , p j2 ) and p ∧ = (p j3 , . . . , p j ℓ ). Similarly, for any vector x ∈ R n , let x ∨ = (x j1 , x j2 ) and x ∧ = (x j3 , . . . , x j ℓ ).

To emphasize that b ∨ and b ∧ depend only on x ∨ and x ∧ , respectively, with an abuse of notation we will write b ∨ (x ∨ ) and b ∧ (x ∧ ) in place of b ∨ (x) and b ∧ (x). For any nonempty subset D ′ of D, define the sets

L max,∨ (D ′ , p ∨ ) = Conv   k∈K(D ′ ) ≤ (x k ∨ , p ∨ )   , (15a)
U max,∨ (D ′ , p ∨ ) = k∈K(D ′ ) ≥ (x k ∨ , p ∨ ) , (15b)
L min,∧ (D ′ , p ∧ ) = k∈K(D ′ ) ≤ (x k ∧ , p ∧ ) , (15c)
U min,∧ (D ′ , p ∧ ) = Conv   k∈K(D ′ ) ≥ (x k ∧ , p ∧ )   . ( 15d
)
The next proposition clarifies the approximation properties of the various sets L and U in (15). ‡ "c-" stands for convex.

Copyright c 0000 John Wiley & Sons, Ltd. Int. J. Robust. Nonlinear Control (0000) Prepared using rncauth.cls DOI: 10.1002/rnc

Proposition 5

Let M ∨ (D ′ ) = max{b ∨ (x k ∨ ) : k ∈ K(D ′ )} and µ ∧ (D ′ ) = min{b ∧ (x k ∧ ) : k ∈ K(D ′ )}. Then, L max,∨ (D ′ , p ∨ ) ⊆ B M∨(D ′ ) (b ∨ ) , (16a) U max,∨ (D ′ , p ∨ ) ⊆ T M∨(D ′ ) (b ∨ ) , (16b) L min,∧ (D ′ , p ∧ ) ⊆ B µ∧(D ′ ) (b ∧ ) , (16c) U min,∧ (D ′ , p ∧ ) ⊆ T µ∧(D ′ ) (b ∧ ) , (16d)
or, equivalently,

x ∨ ∈ L max,∨ (D ′ , p ∨ ) =⇒ b ∨ (x ∨ ) ≤ M ∨ (D ′ ) , (17a) x ∨ ∈ U max,∨ (D ′ , p ∨ ) =⇒ b ∨ (x ∨ ) ≥ M ∨ (D ′ ) , ( 17b
) x ∧ ∈ L min,∧ (D ′ , p ∧ ) =⇒ b ∧ (x ∧ ) ≤ µ ∧ (D ′ ) , ( 17c
) x ∧ ∈ U min,∧ (D ′ , p ∧ ) =⇒ b ∧ (x ∧ ) ≥ µ ∧ (D ′ ) . ( 17d
)
According to [START_REF] De Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF], points in the various sets U and L provide upperbounds and lowerbounds to minima and maxima of b ∨ and b ∧ over D ′ . For example, from (17a) one has that b ∨ (x ∨ ) is a Lowerbound (L) to the maximum (max) M ∨ (D ′ ), hence the notation L max,∨ . Similarly, b ∧ (x ∧ ) is an Upperbound (U) to the minimum (min) µ ∧ (D ′ ), whence U min,∧ . Equivalently, sets U and L provide inner approximations for the various sets T and B. The idea is now to combine these results to establish inequalities for b ∨ (x ∨ )b ∧ (x ∧ ), and hence for the measured values of g

(x) = κ 0 + κ 1 b ∨ (x ∨ )b ∧ (x ∧ ). Recall that µ(D ′ ) = min{g k : k ∈ K(D ′ )} and let M (D ′ ) = max{g k : k ∈ K(D ′ )}. Proposition 6 If g ∈ S 1 (j ∨ , j ∧ , p), for any nonempty set D ′ ⊆ D one has {x ∈ R n : x ∨ ∈ L max,∨ (D ′ , p ∨ ), x ∧ ∈ L min,∧ (D ′ , p ∧ )} ⊆ B M (D ′ ) (g) , {x ∈ R n : x ∨ ∈ U max,∨ (D ′ , p ∨ ), x ∧ ∈ U min,∧ (D ′ , p ∧ )} ⊆ T µ(D ′ ) (g) .
The inner approximations provided by Proposition 6 lead to the following criteria for the invalidation of the family S 1 (j ∨ , j ∧ , p).

Definition 6

The structure (j ∨ , j ∧ , p) is c-inconsistent if there exists a nonempty set D ′ ⊆ D and a data point (x * , g * ) ∈ D \ D ′ such that either of the following conditions applies:

(I) x * ∨ ∈ L max,∨ (D ′ , p ∨ ), x * ∧ ∈ L min,∧ (D ′ , p ∧ ), g * > M (D ′ ) (II) x * ∨ ∈ U max,∨ (D ′ , p ∨ ), x * ∧ ∈ U min,∧ (D ′ , p ∧ ), g * < µ(D ′ ) The structure (j ∨ , j ∧ , p) is c-consistent if it is not c-inconsistent.
Definition 6, that is illustrated in Fig. 3, is a strengthening of Definition 2 to S 1 models. In particular, when |D ′ | = 1, condition (II) is equivalent to the condition in Definition 2. Therefore, if p is m-inconsistent, no g ∈ U (p) (and hence no g ∈ S 1 (p)) can generate the dataset D. However, if p is m-consistent the family S 1 (j ∨ , j ∧ , p) is still invalidated if the structure (j ∨ , j ∧ , p) is c-inconsistent.

Remark 3

Following on Remark 2, if one sets j ∨ = ∅, conditions x * ∨ ∈ L max,∨ (D ′ , p ∨ ) and x * ∨ ∈ U max,∨ (D ′ , p ∨ ) become empty. In this case, condition (II) in Definition 6 coincides with the condition in Definition 5. Moreover, with j ∨ = ∅, it is possible to show that condition (I) in Definition 6 is encompassed by m-inconsistency conditions and hence by condition (II). Therefore, c-inconsistency of the structure (∅, j ∧ , p) is equivalent to c-inconsistency of the sign pattern p.

We highlight that inconsistency conditions in Definition 6 could be easily adapted to address the more general case of models with b(x) = b(x) b(x), where x and x are distinct subvectors of x, b(x) is quasi-concave and b(x) is quasi-convex.

Hierarchical properties of c-inconsistency

Following on Remark 1, for S 0 ∪ S 1 models, it is possible to establish hierarchical relationships among the model structures invalidated by the data. To this purpose, we allow for j ∨ = ∅ and introduce a partial order relating nested model structures. The structure s ′ = (j ′ ∨ , j ′ ∧ , p ′ ) is a substructure of s = (j ∨ , j ∧ , p) (and s is a superstructure of s ′ ) if p ′ is a subpattern of p (see Remark 1), j ′ ∨ ⊆ j ∨ and j ′ ∧ ⊆ j ∧ . Using arguments similar to [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] for the m-inconsistency analysis, it is possible to show that the following properties hold:

• if a structure is c-inconsistent, then all its substructures are c-inconsistent;

• if a structure is c-consistent, then all its superstructures are c-consistent. Such properties have two important consequences. First, they allow one to avoid testing cconsistency of a structure if a substructure (respectively, a superstructure) is already found to be c-consistent (respectively, c-inconsistent). Second, it is possible to provide a compact description of the hierarchy of c-consistent structures by means of its minimal elements (with respect to the substructure partial order). These features are exploited to set up an efficient exploration of all possible structures, as shown in the next section.

ALGORITHMS AND IMPLEMENTATION

An efficient method for testing m-inconsistency was proposed in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF]. The procedure is based on Proposition 1 and hierarchical properties of sign patterns (see Remark 1) rather than the geometric approach discussed in Section 2.3. Here we are concerned with the practical use of Definitions 5-6 for testing c-inconsistency of S 0 ∪ S 1 models. In the following discussion, in view of Remark 3, we will allow j ∨ to be empty and focus on Definition 6 only.

A direct application of Definition 6 is impractical since conditions (I) and (II) must be checked for all subsets of D. In the sequel we show that, without loss of generality, it is possible to check inconsistency of a structure by constructing only two subsets of D for each (x * , g * ) ∈ D.

Proposition 7

A structure (j ∨ , j ∧ , p) is c-inconsistent if and only if there exists (x * , g * ) ∈ D such that either of the following conditions apply:

(I') x * ∨ ∈ L max,∨ (D L , p ∨ )
where

D L = {(x, g) ∈ D \ {(x * , g * )} : g < g * , diag(p ∧ )(x ∧ -x * ∧ ) ≥ 0} ; (18)
Algorithm 1 c-inconsistency test for a structure s = (j ∨ , j ∧ , p) where

D U = {(x, g) ∈ D \ {(x * , g * )} : g > g * , diag(p ∨ )(x ∨ -x * ∨ ) ≤ 0} . (19)
The rightmost inequalities in ( 18) and ( 19) are interpreted componentwise.

The complete method for testing c-inconsistency of a given structure is summarized in Algorithm 1. Note that, when j ∨ = ∅ (respectively, j ∧ = ∅), condition (I') (respectively, condition (II')) in Proposition 7 is not of interest hence one can ignore line 3 (respectively, line 4) of Algorithm 1.

For the efficient implementation of Algorithm 1, it is crucial to have a computationally efficient method for verifying conditions

x * ∨ ∈ L max,∨ (D L , p ∨ ) and x * ∧ ∈ U min,∧ (D U , p ∧ ).
In particular, we would like to avoid computing convex hulls in (15a) and (15d). Inspired by [START_REF] Hanoch | Testing the assumptions of production theory: A nonparametric approach[END_REF], we propose a solution based on Linear Programs (LPs). It is easy to see that conditions (I') and (II') are both instances of the following problem: given points z

1 , . . . , z K , z * ∈ R d and p ∈ {-1, 1} d , check if z * ∈ Z = Conv K k=1 {z ∈ R d : p j z j ≥ p j z k j , ∀j = 1, . . . , d} . (20)
In an equivalent way, ( 20) is false if and only if there exists a hyperplane h T z = h 0 , for some h ∈ R d and h 0 ∈ R, separating z * from the polyhedron Z. It is easily seen that if such a hyperplane exists, then one exists with normal direction aligned with p, i.e. fulfilling p i h i ≥ 0 for all i = 1, . . . , d.

Under this condition, one can just seek a hyperplane passing through z * and such that h T z k > h 0 for all points z k . Testing condition [START_REF] Keller | Model genetic circuits encoding autoregulatory transcription factors[END_REF] then amounts to solve the LP

max δ∈R,h∈R d δ (21) s.t. p i h i ≥ 0, ∀i = 1, . . . , d h T (z k -z * ) -δ ≥ 0, ∀k = 1, . . . , K n i=1 p i h i = 1
and then check if the optimal cost is nonpositive. It is possible to show that ( 21 As already mentioned, if the pattern p is m-inconsistent, all structures (j ∨ , j ∧ , p) are also cinconsistent. However, Algorithm 1 requires the solution of LPs, making c-inconsistency tests more computationally demanding than m-inconsistency tests, for which the efficient algorithm in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] can be applied. Therefore, substantial computational savings can be achieved using the procedure reported in Algorithm 2 for computing the set S min of minimal (with respect to the partial order on structures) c-consistent structures, which exploits the m-inconsistency analysis and the hierarchical properties discussed in Section 3.3.

The correctness of Algorithm 2 (i.e. the fact that the output S min is indeed the set of all minimal consistent structures) can be proven by the same techniques used in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] for the computation of for all structures s = (j ∨ , j ∧ , p) such that p is m-consistent and has ℓ nonzero entries do 5:

if there is no structure s ′ ∈ S min such that s ′ is a substructure of s then 6:

if s is found c-consistent by Algorithm 1 then end for 11: end for minimal consistent sign patterns. In particular, ℓ in line 3 is the number of regulating genes, and hence Algorithm 2 tests c-consistency of simpler structures first. Algorithm 2 was implemented in Matlab 7.10 (R2010a), resorting to the free solver CDD [START_REF] Fukuda | C-library cddlib[END_REF] and its Matlab interface CDDmex [START_REF] Kvasnica | Cddmex -matlab interface for the cdd solver[END_REF] for solving the LP (21).

Handling noisy data

To deal with noisy measurements of (x k , g k ) in D, we follow a robust approach. We assume lower and upper bounds l(•) and u(•) to be available for the true values of g k and x k j , for all k = 1, . . . , m and j = 1, . . . , n. This means that every x k is surrounded by an uncertainty box. The example in the next section shows that this approach is still viable in the case of Gaussian (unbounded) noise affecting the data. The idea is to robustify all inconsistency conditions by defining worst-case scenarios that take bounded uncertainty into account. For what concerns g, conditions g * < µ(D ′ ) and g * > M (D ′ ) are replaced by

u(g * ) < min k∈K(D ′ ) l(g k ) and l(g * ) > max k∈K(D ′ ) u(g k ) , (22)
respectively. Conditions on points x all involve sets computed as combinations (union, intersection, convex hull) of cones. In this case we consider worst-case inner approximations of such sets obtained by replacing ≥ (z, p), ≤ (z, p) with

≥ (z, p) = {z ′ ∈ R d : z ′ j ≥ u(z j ), ∀j such that p j = 1, z ′ j ≤ l(z j ), ∀j such that p j = -1} , (23a)
≤ (z, p) = ≥ (z, -p) = {z ′ ∈ R d : z ′ j ≤ l(z j ), ∀j such that p j = 1, z ′ j ≥ u(z j ), ∀j such that p j = -1} . (23b)
respectively. Moreover, the uncertainty of the test point x * is also taken into account by considering the point x * instead, with

x * j = l(x * j ), j such that p j = 1, u(x * j ), j such that

p j = -1, (24a) or
x * j = u(x * j ), j such that p j = 1, l(x * j ), j such that for conditions involving ≥ or ≤ cones, respectively. As an example, condition x * ∨ ∈ L max,∨ (D ′ , p ∨ ) for p ∨ = (1, 1) is replaced by testing if x * ∨ = (u(x * j1 ), u(x * j2 )) belongs to the set

p j = -1, (24b)
Lmax,∨ (D ′ , p ∨ ) = Conv k∈K(D ′ ) ≤ (x k ∨ , p ∨ ) , (25)
as represented in Fig. 4, to be compared to Fig. 3(a). We highlight that, for conditions involving convex hulls, this approach is equivalent to replacing the LP [START_REF] Sontag | Network reconstruction based on steady-state data[END_REF] with

max δ∈R,h,z 1 ,...,z K ,z * ∈R d δ s.t. p i h i ≥ 0, ∀i = 1, . . . , d h T (z k -z * ) -δ ≥ 0, ∀k = 1, . . . , K n i=1 p i h i = 1 l(z * ) ≤ z * ≤ u(z * ) l(z k ) ≤ z k ≤ u(z k ), ∀k = 1, . . . , K
Roughly speaking, larger uncertainties, that correspond to higher u(x k j )l(x k j ) and u(g k )l(g k ), result in a lower chance that a structure is declared c-inconsistent. In particular, ( 23)- [START_REF] Fukuda | C-library cddlib[END_REF] shrink sets L and U in [START_REF] Kauffman | Genetic networks with canalyzing boolean rules are always stable[END_REF] and shift the test point x * , resulting in a smaller number of test points that can be used for falsifying a structure.

EXAMPLES AND DISCUSSION OF THE RESULTS

In order to assess the falsification capability provided by quasi-convexity properties, we considered the same artificial network introduced in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] for evaluating the performance of the m-inconsistency analysis. The network, represented in Fig. 5, comprises 6 genes and several interactions. In particular, genes 1-3 represent the core oscillating part of the system and correspond to the repressilator network developed and synthesized in Escherichia coli [START_REF] Elowitz | A synthetic oscillatory network of transcriptional regulators[END_REF]. The remaining three genes are those of interest in our study and are regulated by the three core genes according to different logical rules. The dynamics of this part of the network is modeled by

ẋ4 = κ 0,4 + κ 1,4 σ -(x 1 )σ + (x 2 ) -γ 4 x 4 , (26a)
ẋ5 = κ 0,5 + κ 1,5 [1 -σ + (x 2 )σ -(x 3 )] -γ 5 x 5 , (26b)
ẋ6 = κ 0,6 + κ 1,6 [1 -σ + (x 2 )σ + (x 3 )]σ + (x 1 ) -γ 6 x 6 . (26c)
where x i , i = 1, . . . , 6, denotes the concentration of the ith gene product and σ ± (•) are Hill functions. We defer the reader to the supplementary material of [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] for the details about the complete model and the parameter values. It is easy to recognize regulation functions in both S 0 (gene 4) and S 1 (genes 5 and 6).

We have tested the performance of both the m-inconsistency analysis introduced in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] and its combined use with c-inconsistency analysis developed in this work on the following datasets. The model was simulated for 15 time units § , i.e. until 3 full oscillations were completed. To evaluate the sensitivity of our approach to the amount of data, we produced three datasets comprising m = 45, 23 and 12 equally spaced data points. In order to assess the impact of measurement errors, noisy synthesis rate and concentration data gk and xk were obtained by corrupting g k and x k with multiplicative noise, i.e., for i = 1, . . . , n,

xk i = x k i (1 + s e e k i ) , (27)
gk i = g k i (1 + s ǫ ǫ k i ) , (28)
with e k i and ǫ k i mutually uncorrelated Gaussian random variables with zero mean and unit variance [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF]. The scaling factors s e , s ǫ were chosen in the set {0.03,0.05,0.07}. Recalling that, approximately, noise samples s e e k i fall within ±3s e with 0.99 probability (and similarly for s ǫ ǫ k i ), we considered noise contributions ranging from the 9% to the 21% of the noiseless data values. In order to establish lower and upper bounds to the data, as assumed in Section 4.1, we used 95% confidence intervals resulting in

l(x k i ) = xk i (1 -2s e ) , u(x k i ) = xk i (1 + 2s e ) , (29)
l(g k i ) = gk i (1 -2s ǫ ) , u(g k i ) = gk i (1 + 2s ǫ ) . (30)
Since this gene network model involves Hill functions, the log-transformation ( 7) was applied to the data prior to the execution of Algorithm 1. We are interested in comparing the falsification performance of the c-consistency analysis to the case when only m-inconsistency is used, i.e. when only the method of [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF] is applied. To this purpose, let N c-inc be the number of c-inconsistent structures and N m-inc the number of minconsistent structures. The performance index I % ≥ 0 is defined as and the larger the I % , the larger the percentage of S 0 ∪ S 1 models invalidated by c-inconsistency but not by m-inconsistency that is, the larger the increase in performance with respect to the method in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF]. In order to quantify the fraction of all structures that are falsified, we also introduce the selectivity index

S % = N c-inc |S 0 ∪ S 1 | • 100 , (32)
where the total number of structures is |S 0 ∪ S 1 | = 5588 for the considered network. The larger S % the larger the portion of S 0 ∪ S 1 structures invalidated by c-inconsistency. Average values of the performance indices are reported in Table I for varying values of dataset size m and noise scaling factors s e , s ǫ . They were obtained from 100 Monte Carlo experiments, each characterized by different noise realizations. The true structure was never declared inconsistent, showing the robustness of the falsification procedure. This suggests that, despite the fact that bounds l(•), u(•) do not take into account 5% of variability, the falsification conditions are rather robust for the dataset considered in our study. Concerning the selectivity index S % , one can notice a degradation of performance when either the noise level increases or the size of the dataset decreases. In both cases, this is due to the fact that less datapoints can be used for model falsification. While obvious for smaller datasets, with the increase of noise this behavior is explained as follows. Refer for instance to Fig. 4. Model invalidation relies on checking whether x * ∨ belongs to set Lmax,∨ . For larger values of noise, bounds are increased so that x * ∨ is robustly classified (with high probability) inside or outside L max,∨ based on the noisy observations. This conservatism, which ensures that the true model is not invalidated (with high probability), comes at a price: less datapoints are appropriate for invalidation, whence the decrease of S % . The variability of S % among the three genes also suggests that the considered datasets do not equally support structure falsification for different genes. The analysis of index I % highlights an interesting behavior. Excluding the least favorable condition s e = s ǫ = 0.07 and m = 12, the contribution of the c-inconsistency analysis increases when datasets become smaller and noisier. A more detailed analysis of the results reveals that the improvement is most significant (up to 40%) for candidate models with largest number of regulating genes. This means that c-inconsistency can play a key role when structures of high complexity need to be falsified on the basis of few noisy data.

CONCLUSIONS

In this paper, we introduced and analyzed geometrical properties of a relevant class of gene network models. Under the assumption that measurements of gene product concentrations and synthesis rates are available, we exploited monotonicity and convexity-like properties to invalidate families of models that are inconsistent with the data. The proposed falsification techniques represent an extension of the approach presented in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF], where the concept of sign pattern was introduced to capture the monotone character of unate regulation functions. We highlight that the approach proposed in this work can be extended to deal with any regulation function in the form of a product of a quasi-concave and a quasi-convex function. Since this extension would allow for many more alternative model structures, an increase in the overall complexity of the falsification strategy is expected.

The performance of the method was evaluated by means of Monte Carlo simulation using an oscillating synthetic network model. The results demonstrate a substantial improvement with respect to the approach in [START_REF] Porreca | Identification of genetic network dynamics with unate structure[END_REF], especially when a small, noisy dataset is used.

Future directions of this research include the application of the proposed method to real experimental data to confirm the results obtained in silico. On the theoretical side, we anticipate a detailed study of other model classes that can benefit from our invalidation methods as well as the development of criteria for the invalidation of more general subclasses of unate models.

Figure 1 .

 1 Figure 1. Cones ≥ (x 1 , p) and ≥ (x 2 , p) composing the set M (D ′ ) for p = (-1, -1), D ′ = {(x 1 , g 1 ), (x 2 , g 2 )} and g 2 = µ(D ′ ). Dashed lines show the boundaries of super-level sets T g 1 (g) and T g 2 (g).

Figure 2 .

 2 Figure 2. Sets C(D ′ , p) (gray area) and M (D ′ , p) (subset of gray area left of the dotted lines) for p = (-1, -1), D ′ = {(x i , g i ), i = 1, . . . , 4} and g 4 = µ(D ′ ). The dashed line shows the boundary of the superlevel set T g 4 (g). Point x # represents a test point that allows both m-and c-inconsistency falsification, whereas x * allows for c-inconsistency falsification only.

  Set Lmax,∨(D ′ , p∨) in the x∨-plane. (b) Set L min,∧ (D ′ , p∧) in the x∧-plane.

Figure 3 .

 3 Figure 3. Sets in condition (I) of Definition 6 for x ∈ R 4 , D ′ = {(x 1 , g 1 ), (x 2 , g 2 )} and p∨ = p∧ = (1, 1).Dashed lines indicate the boundaries of level sets of b∨ (left) and b∧ (right).

  ) is always feasible, with the last constraint ensuring boundedness. By the use of the LP (21), one execution of Algorithm 1 amounts in the worst case to solving |D| LPs in |j ∨ | + 1 variables and |D| LPs in |j ∧ | + 1 variables, with one equality constraint and at most |j ∨ | + |D| -1 and |j ∧ | + |D| -1 inequality constraints, respectively.

Algorithm 2

 2 Computation of the set S min of minimal c-consistent structures 1: initialize S min = ∅; 2: assess m-inconsistency of sign patterns by means of the algorithm in [1]; 3: for ℓ = 1, . . . , n do 4:

7 :S

 7 min = {s} ∪ S min ;

Figure 4 .

 4 Figure 4. Uncertainty boxes (dashed), set Lmax,∨(D ′ , p∨) in (25) and robustified test point x * ∨ for D ′ = {x 1 , x 2 } and p∨ = (1, 1).

Figure 5 .

 5 Figure 5. Artificial regulatory network: repressilator loop (genes 1-3) plus controlled genes 4-6.

  D L as in[START_REF] Plahte | A methodological basis for description and analysis of systems with complex switchlike interactions[END_REF]. If x * ∨ ∈ L max,∨ (D L , p ∨ ) label s as c-inconsistent and exit. D U as in[START_REF] Yang | An analytical rate expression for the kinetics of gene transcription mediated by dimeric transcription factors[END_REF]. If x * ∧ ∈ U min,∧ (D U , p ∧ ) label s as c-inconsistent and exit.

	5: end for
	(II') x * ∧ ∈ U min,∧ (D

1: label s as c-consistent 2: for all (x * , g * ) ∈ D do 3: compute 4: compute U , p ∧ )

Table I .

 I Performance results on the example network.

			se,sǫ	se,sǫ	se,sǫ
	m gene	0.03	0.05	0.07
		4	S % =47.60 S % =42.30 S % =37.06 I % = 5.07 I % = 5.35 I % = 6.47
	45	5	S % =41.46 S % =31.23 S % =23.84 I % = 6.65 I % = 8.54 I % =12.08
		6	S % =40.19 S % =36.36 S % =31.25 I % = 6.01 I % = 8.30 I % =10.53
		4	S % =45.00 S % =38.57 S % =32.86 I % = 7.92 I % = 8.73 I % = 9.28
	23	5	S % =34.81 S % =26.39 S % =19.03 I % = 8.65 I % =12.13 I % =13.03
		6	S % =37.66 S % =32.86 S % =27.29 I % =10.09 I % =12.97 I % =14.18
		4	S % =41.78 S % =33.37 S % =25.28 I % = 9.63 I % =10.84 I % =10.99
	12	5	S % =29.44 S % =22.37 S % =16.62 I % =11.54 I % =12.76 I % =12.19
		6	S % =32.82 S % =27.88 S % =23.44 I % =10.88 I % =12.87 I % =12.45

Copyright c 0000 John Wiley & Sons, Ltd.Int. J. Robust. Nonlinear Control (0000) Prepared using rncauth.cls DOI: 10.1002/rnc

† "m-" stands for monotone.Copyright c 0000 John Wiley & Sons, Ltd.Int. J. Robust. Nonlinear Control (0000) Prepared using rncauth.cls DOI: 10.1002/rnc

ACKNOWLEDGEMENT

This work was partially supported by the European Commission under the Network of Excellence HYCON2, contract number FP7-ICT-257462, and by the SystemsX.ch research consortium under the project YeastX.


J Bender

A Kuijper

D W Fellner

É Guérin

Marie Durand

Bruno Raffin

François Faure

Workshop on Virtual Reality Interaction and Physical Simulation VRIPHYS (2012)

Neighbor identification is the most computationally intensive step in particle based simulations. To contain its cost, a common approach consists in using a regular grid to sort particles according to the cell they belong to. Then, neighbor search only needs to test the particles contained in a constant number of cells. During the simulation, a usually small amount of particles are moving between consecutive steps. Taking into account this temporal coherency to save on the maintenance cost of the acceleration data structure is difficult as it usually triggers costly dynamics memory allocations or data moves. In this paper we propose to rely on a Packed Memory Array (PMA) to efficiently keep particles sorted according to their cell index. The PMA maintains gaps in the particle array that enable to keep particle sorted with O(log 2 (n)) amortized data moves. We further improve the original PMA data structure to support efficient batch data moves. Experiments show that the PMA can outperform a compact sorted array for up to 50% element moves.

Introduction

Particle systems are commonly used for simulating fluids, crowds, molecular dynamics, granular materials as well as impacts on concrete structures. These methods rely on particles whose state is updated according to interactions with neighbor particles (we consider here only short range interactions). Neighbor search is a computationally intensive step. Acceleration data structures are classically used to reduce its cost. A common approach consists in partitioning the simulated space with a regular grid, identifying the particles each cell contains. Cell size is defined according to the interaction range to limit neighbor search to the particles contained in the local and neighbor cells.

Given the large number of particles a simulation may require, performance can be significantly impacted by the memory layout. Having particles belonging to the same cell close by in memory enables to preserve a low cache miss ratio. But maintaining this locality as particles move may require costly memory copies or dynamic allocations. The challenge is thus to find a data structure that offers a good tradeoff between particle locality and maintenance cost. To reach this goal it is important to take into consider-ation the strong temporal locality that usually show particle simulations: only a small ratio of particles actually change of cell between consecutive iterations.

In this paper we propose to store particles in a Packed Memory Array (PMA). Particles are stored sorted (according to their cell index) and evenly distributed in an array larger than the actual number of particles. To further improve spatial coherency, cells are indexed following a Z-curve pattern. An amortized scheme enables to maintain the PMA in a coherent state at a low cost. Inserting a new particle while maintaining the PMA state has a O(log 2 (n)) amortized cost. Range query or PMA scanning leads to a controlled cache miss overhead proportional to the PMA density. To optimize performance we adapt the PMA to move particles by batches, further reducing the cost of the PMA maintenance. Experiments show that the PMA can outperform a compact sorted array up to 50% elements moving amongst 10M ones. Performance is improved by a factor of 9 for 1% moves and by 50% for 50% moves. Experiments with a Smoothed Particle Hydrodynamics (SPH) simulation confirm the interest of the PMA data structure for keeping moving particles sorted.

The paper is organized in 6 sections. After reviewing related work (Sec. 2), we remind some basics about particle simulation and give an overview the PMA integration (Sec. 3). The PMA data structure is detailed (Sec. 4) with experimental results on random integers. An SPH simulation is then used to probe the PMA for particle simulation (Sec. 5). The conclusion summarizes our contribution and discusses future research directions (Sec. 6).

Motivation and Related Work

The naive approach for computing particle neighbors tests each particle pair, leading to O(N 2 ) tests. The Verlet list [START_REF] Fomin | Consideration of data load time on modern processors for the verlet table and linked-cell algorithms[END_REF] improves this approach by using a neighborhood distance larger than the actual cut-off distance required by the simulation. Thus, the Verlet list, storing the neighbor list of each particle, does not need to be updated at every time step, amortizing its maintenance cost. Using a uniform grid to partition the simulation space, particles are distributed in each cell. Filling the grid is O(N ). Finding a particle's neighbors is O(1) assuming that each cell can only contain a constant maximum number k of particles [START_REF] Fomin | Consideration of data load time on modern processors for the verlet table and linked-cell algorithms[END_REF]. To avoid costly memory allocations, usually an array of size k is allocated for each cell, making it a O(k * m) size data structure, m being the number of cells. All cells, even empty ones, are stored, impairing memory consumption for very large simulation domains [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF].

Spatial Hashing proposes to improve memory consumption [THM * 03]. Particles are distributed in a hash map according to their cell index. Filling the hash map is thus O(N ) and finding a particle's neighbors is O(1). The authors note that having 2 * N entries in the hash function enables to have a low collision rate. To avoid numerous dynamics memory allocations, the non empty cells are stored in a compact array with room for storing k particles per cell [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF]. Index sorting consists in sorting the particles in an array according to their cell index [START_REF] Müller | Particle-based fluid simulation for interactive applications[END_REF]. The construction cost is thus O(N ). Memory usage is limited to the actual number of particles. If all cells store a pointer, set to null for empty cells or pointing to the cell first particle, finding a particle's neighbors is O(1), but simulating large domains can be memory consuming. Memory can be traded for computations by storing only non empty cells, and performing a O(log(N )) binary search to find the neighbor cells of a given particle.

Given the growing memory complexity of modern processors, being CPUs or GPUs, data layout in memory can significantly affect the cache hit or coalesced access ratios [START_REF] Tchiboukdjian | Binary Mesh Partitioning for Cache-Efficient Visualization[END_REF]. Sorting particles to keep spatial locality, i.e. keeping the spatial locality when mapping the data in memory, can significantly decrease the required number of memory transfers. The Z-curve, also called Morton order, or Hilbert order are two classical memory layout patterns that show good locality properties. The Hilbert order often outperforms the Z-curve in term of locality preservation, but at the price of a more complex index computation. On the contrary, the Z-index can be efficiently computed by bitwise operations on the 3D coordinates [START_REF] Butz | Alternative algorithm for hilbert's space-filling curve[END_REF]. For particles simulations, the Z-curve is often used to improve the probability that data related to neighbor cells are stored close by and then are likely to be accessed consecutively. Thus, index sorting classically rely on a Z-curve. Ihmsen et al. [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF] improve spatial hashing, called compact hashing, by storing the actual particle data in a Z-order sorted array. In both cases, sorting is usually not performed at each iteration to amortize the sorting cost. Experiments in [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF] show performance improves significantly when relying on a Z-curve.

Particle simulations usually show a high temporal locality, i.e. only a small ratio of particles move between cells in consecutive time steps. Ihmsen et al. [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF] mention a 2% move ratio, while we experienced up to 25% moves with the fluid simulator [START_REF] Hoetzlein | Fluids v2.0, open source, fluid simulator[END_REF]. One idea is to take advantage of this locality to build specific data structures that support an efficient incremental update, while still ensuring a fast neighborhood search and a low cache miss ratio.

The basic approach often reported is to rely on classical data structures and simply sort the particles periodically. For instance in [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF], particles are sorted once every 100 iterations, leading to a good trade-off between sorting cost and cache efficiency drop as the simulation shows only a 2% move ratio. As our experiments will show (Fig 6 ), sorting at every time step can be beneficial for higher move ratios.

[AE97] presents a general approach to make range query data structures dynamics. The idea is to handle several structures of varying 2 i sizes. When a new element needs to be inserted, all elements gathered in the smallest full j data structures are gathered in a 2 (j+1) size one. The cost of data structure building is thus amortized: the more elements the less often the data structure is rebuilt. This scheme is adapted to insertions but does not support efficient deletions. Following the same principle and inserting gaps amongst the elements of a sorted array, Bender introduced the library sort [START_REF] Bender | Insertion sort is o(n log n)[END_REF]. To support efficient element removals as well as insertions in a sorted array, the data structure was further modified leading to Packed Memory Arrays (PMA) [START_REF] Bender | Cache-oblivious b-trees[END_REF][START_REF] Bender | An adaptive packedmemory array[END_REF]. A controlled gap density is maintained into an array of sorted elc The Eurographics Association 2012.

ements enabling to insert or delete an element in O(log 2 (N )) amortized element moves. These gaps obviously require extra storage space and generate additional memory transfers during scans. The PMA thus enables to control the trade-off between maintenance cost and memory dispersion. PMAs have been developed for large databases in conjonction with Cache Oblivious B-Trees [START_REF] Bender | Cache-oblivious b-trees[END_REF], and, to our knowledge, never used for numerical simulations. We propose in this paper to keep particles Z-sorted using a PMA. In simulation, contrary with databases, the updates are applied in large batches at the end of each time step. Thus we extended the PMA to support batch insertions, deletions and moves.

Overview

Before detailing the PMA and its application to particle sorting, we review some of the basics of particle simulations and give an overview of the PMA use with particles. Particle simulations classically perform first a neighbor search, based on a cut-off distance. The entire list of particles is next scanned several times to update particle states. For most particle-based simulations, (Alg. 1), densities and next forces are updated based on particle neighborhoods, before the integration and position update that only require to access particle local data. A naive neighborhood search is thus O(N 2 ) while all other scan steps are O(N ).

Algorithm 1 Particle simulation update State neighborhood search for all particles do compute density, pressure from neighborhood for all particles do compute force from interaction with neighbors for all particles do update positions

As mentioned in previous section, a classical approach to accelerate neighborhood search is to divide the simulation space into a uniform 3D grid of cells. Each particle is assigned to the cell its centroid belongs to. Finding the neighbor particles only requires to test the particles belonging to the current and neighbor cells.

Particles belonging to the same cell can simply be linked. Though it avoids data movements, it can seriously affect cache performance. The particles of a cell are accessed together. If the particles of a cell are spread in memory, accessing them may potentially lead to load numerous cache lines with a low relevant content (one particle per cache line in the worst case). To avoid this issue, particles can be sorted according to their cell index (Fig 1). To further improve memory locality, cells can be indexed following a Z-curve. Thus the particles of spatially close cells will tend to be stored nearby in memory. Though cache efficiency is increased, sorting particles introduces an extra cost. With an array of contiguous particles, moving one of them can lead to shift every other element, whatever the method we choose (re-sorting the whole array with standard sorting algorithms, or moving particles using an insertion sort like pattern).

Inside an array, the only way to reduce the amount of memory moves is to anticipate insertions keeping extra (empty) spaces. With enough gaps, there is always an empty space close to the target insertion location that can be used to make room for the element to insert with a low number of memory moves. The space density should nevertheless be limited to avoid increasing the memory consumption and memory transfers. Indeed, these gaps will obviously make the overall array bigger and decrease the ratio of the relevant data a loaded cache line can contain.

The PMA structure implements this apparently simple idea with an amortized scheme for keeping a bounded gap density as elements are removed or inserted. As already mentioned we propose to use the PMA to store particles and keep them sorted according to their cell index.

PMA

We present the PMA algorithm as introduced in Θ(N/log N ) pieces called segments (Fig. 2). The density of segments, i.e. the ratio of elements against gaps, is maintained between the density bounds ρ0 and τ0.

The density of a pair of consecutive segments, level l = 1, is maintained between ρ1 and τ1. This scheme repeats by doubling the number of segments considered up to include the full array bounded by densities ρ h and τ h . Densities are linearly staged between the different levels. At level l, the density of a window of 2 l segments is bounded by:

τ l = τ h + (τ0 -τ h ) (h -l) h (1)
ρ l = ρ h + (ρ0 -ρ h ) (h -l) h (2)
The density gap between two consecutive levels l and l + 1 is O(1/ log N ) to ensure the complexity of the algorithm [START_REF] Bender | An adaptive packedmemory array[END_REF].

Extra constraints on maximum and minimum densities ensure that:

ρ0 < ... < ρ h < τ h < ... < τ0 (3)
Thus, if a window of size l respects its density bounds, then all sub-windows also respect their density bounds.

When inserting a new element inside a PMA (Alg. 2), the destination segment is first found. Before proceeding to the insertion, we check that the target segment density, counting the element to insert, respects the segment upper density bound τ0. If not, we look for the smallest windows at level l enclosing the segment that respects its upper density bound τ l . We next proceed to an operation called rebalance that shares out the gaps among the different segments in that window and insert the new element. A rebalance consists in two steps. First elements are packed at one end of the window. Then they are moved into their destination segment, by one scan of the concerned segments. The obtained PMA respects all density bounds.

If no window can be found, i.e. τ h is not respected, the PMA size is doubled and rebalanced. Element deletion is similar except that the ρ l densities must be considered. If the whole PMA becomes too empty, i.e. ρ h is not respected, the PMA size is divided by two. To make sure that when halving or doubling the array size, the new PMA is within the density bounds, the extra constraint on densities must be respected: The PMA is a cache-oblivious data structure: the PMA is built without taking into account a given cache line or total cache size. The performance is anyway guaranteed, with O(1 + S/B) memory transfers when scanning S consecutive elements on a machine with a B size cache line. Cache-oblivious data structures have the advantage, over cache-aware ones, to be more flexible. For instance, they usually can be efficiently shared between a GPU and a CPU.

2ρ h < τ h (

Batch Moves

We now present how we extended the PMA to improve the performance of moving batches of elements. An element is moving when its key, i.e. cell index, is changed. A move can be implemented by first removing and next inserting this element. However this may lead to rebalance twice the same window inside the PMA. When several elements are moving at the same time, the chances that some array section get rebalanced several times increases. To save memory movements we propose to apply moves batchwise.

First, moving elements are removed from the PMA without rebalance. A gap is simply set where the element to remove is. Moving elements are copied in c The Eurographics Association 2012. a temporary array, called moving array. All moving elements are then processed together.

The moving array is sorted using a classic sort algorithm like qsort. We then consider the first element of the PMA middle segment. The position of this pivot is searched (binary search) in the moving array. The moving elements on the left of the pivot need to be inserted in the first half of the PMA, the others in the right part. We then test the density against its allowed bounds, counting the elements to insert, for each of these 2 windows delimited by the pivot.

If the window density does not respect the bounds, its rebalance is triggered. The elements to insert as well as the element in the PMA are already sorted. The rebalance can then be implemented in a single scan on the window, using an extra temporary array to store elements that need to be shifted forward to leave room either for an element to insert or to a gap.

If the density bounds are respected, the process is recursively applied down to individual segments as long as no rebalance is triggered in a larger window. As long as elements need to be inserted, we need to check the window density against its upper and lower bound. Since moving element were previously removed from the PMA without rebalance, some window densities may be below the lower threshold, even if some elements need to be inserted. Once no more element need to be inserted in a window, only the lower density bound needs to be probed. This approach equally applies to raw suppressions or new insertions. New elements to insert are sorted with the moving elements in the moving array. The recursive process starts one level higher, testing the density of the full PMA to check if halving or doubling the PMA size is necessary.

In case of random uniform moves, rebalance will take place on the smallest windows, i.e. at segment level, with a high probability. The move cost is dominated by the sorting of the moving array and the rebalance of each segment. For a high move ratio, it is probably more efficient to simply sort all element in a dense array and next spread them in a PMA, thus saving the overhead associated with the recursive descent along the PMA windows.

Experimental Results

Density Bounds and Segment Size

The size of the PMA, N , is the smallest power of 2 greater than or equal to K/τ h , where K is the number of valid elements to be sorted in PMA. The number of segments, S, is the smallest power of two larger than N/ log(N ) as suggested in [START_REF] Bender | Insertion sort is o(n log n)[END_REF] of a segment is given by N/S. The minima densities for the segments (ρ0) and the whole array (ρ h ) are chosen by symmetry around 0.5 from the maxima ones (τ0 and τ h resp.). The densities ρ h and τ h must also comply with the constraint 2ρ h < τ h (Eq. 4).

No explicit parameter space exploration was performed but we tested various maxima densities τ0 and τ h for various K values. We observed in particular that given a number K of elements, τ h should be chosen such as the final density is between 50% and 95%. When the final density is lower than 50%, there are too many gaps with respect to the number of elements, elements are then spread out accross a larger memory space, creating more memory transfers. Moreover, the size of the PMA being bigger (actually, the double) than the optimal one, it leads to more segments, thus extra PMA maintenance. On the other hand, if the final density is higher than 95%, the segments are almost all full. The insertion is more likely to cause rebalance at a high level.

In the experiences presented in next sections, we take τ0 = 0.92, τ h = 0.7, ρ h = 0.3, ρ0 = 0.08 as these values gave good results for different PMA sizes.

PMA Evaluation

We compare the performance of the PMA against a dense array sorted using a quicksort and an insertion sort for different percentages of moves. All data structures are initialized with random sorted integers (uniform distribution). A defined percentage of elements are randomly elected and assigned new (random) keys. The dense array is sorted again with a quicksort or insertion sort. A batch move is executed on the PMA. We measure the time taken to update all these data structures.

The CPU used is a 4-core Intel Xeon E5520 @2.27GHz (only one core is used for running a simulation) with 18Go of memory. The compiler used is g++-4.6.2. We average execution times over 32 passes.

We use the highly optimized implementation of the quicksort algorithm provided by the libc, qsort, and a home made insertion sort code, isort. Figure 3: Speed-up of PMA versus qsort against % of moves, for 100 000, 1 000 000 and 10 000 000 elements. The speed-up is the ratio of the sort execution time using qsort to the sort execution time using the PMA.

and the insertion sort. The insertion sort, often recognized as the most efficient sort for almost sorted arrays, is not competitive even with the quick sort. The insertion is efficient for almost sorted arrays for K small. But for large K values as the one tested here, the cost of data movements becomes prohibitive. Thus experiments were not performed for insertion sort with K ≥ 1 000 000.

Figure 3 presents the results for a move ratio ranging from 1% to 50% on 100 000, 1 000 000 and 10 000 000 integers. The results show that with 10 000 000 elements, the PMA is about 9 times faster than qsort for 1% moves. With 50% moves, the speedup is 2. As expected, the PMA performance decreases along with the percentage of moves. Also notice that the PMA efficiency increases along with the number of elements.

PMA Overhead

If gaps enable to efficiently maintain moving elements sorted, they also introduce an overhead during range queries. Compared to a dense array, the valid elements are spread on a larger memory region, causing extra memory transfers, and a check is required at every location to probe the element validity (validity check ). We evaluated this overhead through a scan operation on all the elements of a PMA in the following way: we build a PMA with a set of random elements and we benchmark the time to compute the sum of all the elements. We compare to the time needed for summing the same elements stored in sorted dense array. A direct implementation (one loop iterating on the PMA and a validity check for each element) leads to a poor performance. Beyond the expected overhead, we suspect that the validity check impairs some compiler optimizations and cause stalling in the processor pipe-line. After experimenting several approaches, the most efficient one consists in storing the element validity in a bitmask array. When reading one element of this bitmask array, we get 64 bits encoding the validity of 64 PMA elements. The scan is then written using 2 loops, an outer loop iterating on the bitmask elements, i.e. considering 64 PMA elements per iteration, and an inner loop probing the validity of these 64 elements. For K = 1 000 000 elements (Table 2), stored in a PMA array of size N = 2 097 152 = 2.09K, the scan is 3.7 times slower than the scan on a dense array. For K = 10 000 000 elements, the size of the PMA is N = 16 777 216 = 1.7K, explaining the improved performance ratio (3.1 for a fill ratio of 59% instead of 3.7 for 48%).

K array

The sum operation performed per element is lightweight. More time consuming operations will make the overhead less noticeable. If many range queries operations need to be performed, like for a fluid simulation, other approaches can be considered to amortize this overhead. As we discuss in Section 5.2, we build an indirection array pointing to the PMA valid elements to save on the validity check.

Particle Simulation with the PMA

Next step is to use the PMA for a particle simulation. We tested the PMA in Fluids [START_REF] Hoetzlein | Fluids v2.0, open source, fluid simulator[END_REF], an Open Source fluid simulator. We detail implementation and performance issues in the following.

The native Fluids implementation relies on a double array, one storing all cells, even empty ones, and a second one storing particles. By default particles are not sorted. Fluids fills the cells through a scan, each cell pointing to its first particle, while the remaining particles of the cell are chain linked, this operation is called particles insertion in grid.

We replaced the particle data array by a PMA that takes care of maintaining the particles sorted. To avoid c The Eurographics Association 2012. the creation of a big array for all cells of the simulation space, we can also use a PMA to store non empty cells (Fig. 4). The batch insertions and deletes described in section 4.1 are used to keep up to date the non empty cells. The sizes of the arrays are thus only dependent on the particle number, and not on the simulation domain size.

We implemented a Z-curve cell indexing used to sort the cells and particles. A parameter enables to control the frequency of particle sorting for the native Fluids code, while the PMA maintains the particles sorted at each iteration.

As mentioned by [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF], sorting particle indices first and then updating accordingly the actual particle data can improve performance. This optimization could be implemented for both cases, but the experiments presented here rely on direct sort of particle data.

Fluids originally uses a cell size of 2 times the influence radius h. This choice allows to test only 8 neighbor cells per particle, provided that we take into account the position of the particle inside the cell.

We implemented a more classical grid, with a cell size of h that requires to test for proximity the particles of the 27 neighbor cells. We can further reduce the number of particle pairs to test by taking into account the symmetry of interactions. Indeed, we can test only 14 neighbor cells to find only once each pair. In that case, we need to update the neighbors of both interacting particles. Table 3 presents experimental results comparing the different methods for a scene of 130 000 particles, using the standard contiguous arrays CD from Fluids. The collision detection execution time can be divided by 2 when testing only the particles of 14 cells compared to the initial Fluid implementation testing the content of 8 larger cells.

Performance Results

Figure 5 presents the execution time of the sort operation during a corner breaking dam (CBD) simulation with 2 900 000 particles and 1 830 180 cells (354x94x55). The PMA performance is compared to the fluid simulation code in the situation where particles are sorted at every time step, called reference implementation. All cells (empty or not) are stored in a dense array.

The reference sort implementation is organized in four steps. A first pass builds an array of pairs consisting of every particle position in memory and their key according to the Z-curve. Qsort is then applied to that array and not on the whole data array. Then a temporary data array is used to copy the particles to their target position in memory. Finally, we copy the temporary array to the main particle data array. This saves time because the sort itself is, in this case, applied to smaller element sizes. In spite of this optimiza- tion, which benefits only to the reference implementation, the PMA is three times faster than the reference implementation for about 10% of particle moves. Performances are similar around 20% of particle moves.

The corner breaking dam simulation starts with particles falling down till they reach the ground. At this time, particles move quickly, leading to a high move ratio. With the stabilization of the system (here after about 100 frames), the number of moves decreases. We globally experienced with the Fluid simulator a high level of particle moves compared to the 2% discussed in [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF]. Particle simulations with such a low level of particle move will significantly benefit from using a PMA data structure.

PMA Overhead Impact on the Whole Simulation

Table 4 presents the relative execution times of the different steps of a CBD simulation of 180 000 particles.

Here the collision detection execution time is an order of magnitude higher than the sort execution time. This is again a noticeable difference with [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF], where the authors even discuss the possibility of recomputing neighbors several time per time step rather than storing them. Such approach is of course not relevant with such a high cost for collision detection.

This high collision detection cost versus sorting time is also not favorable to the PMA. The overheads absorb all the benefit of the fast particle sorting. The timing of the different steps reveals that PMA overhead is almost not noticeable for the pressure and force computations neither for position updates. The overhead is small compared to the time consuming operations executed for every particle. However, we initially experienced a high overhead during the collision detection step. To reduce this cost, we build an indirection array of PMA valid positions before the beginning of the neighborhood search. The goal is to remove the validity test that can prevent compiler optimizations and filling the processor pipe-line, and replace it with an access to an indirection array pointing to the valid elements in the PMA.

If particles are sorted in memory, it is easy to insert particles in the different cells in one single scan because particles of the same cell are contiguous in memory. Thus we can use the indirection array to store only the first and the last particle of each non empty cell. This latter improvement is suitable as well for Fluids as long as the particles are sorted in memory. It saves 12% of time for the whole collision detection (particles insertion in grid and neighbor search) for Fluids. Combined with the use of the indirection array, it saves 35% for the PMA.

Figure 6 presents the total execution time for simulation with 180 000 particles and 132963 cells (141x41x23). The PMA performance is compared to the native fluid simulation when we never sort particles in memory (no sort), when particles are sorted every 100 time steps and when they are sorted every time step. It is clear that sorting particles improves significantly the performance of the native Fluids code, showing the benefit of maintaining a good memory locality. On the other hand, for the whole simulation, sorting every time steps saves 5% of the computation time against sorting every 100 time steps. This is another difference we noticed compared to [START_REF] Ihmsen | A parallel sph implementation on multi-core cpus[END_REF] where they get the best trade-off between cache efficiency and sort overhead by sorting only every 100 time steps. Finally, using the PMA gives the best performance, saving 4 extra % of the total execution time.

Conclusion

In this paper, we proposed to store particles in a Packed Memory Array. A PMA maintain gaps be- : Total execution times (s) for a CBD simulation of 180 000 particles with native code when particles are never sorted in memory after initialization (No Sort), when they are sorted every 100 time step (Sort Every 100), every time step (Reference), and using the PMA. On the left, we show the execution time for 1 000 time steps. On the right, we focus on the 580 -720 range to highlight differences between sort politics.

tween particles, reducing significantly the cost of moving particles. The PMA significantly outperforms a classical quick sort for a reduced ratio of particle moves, a common situation for particle simulations that show a high level of temporal coherency. The PMA thus appears as a trade-off between a dense array storage leading to costly element moves, and linked arrays requiring either a large, often unused, amount of memory, or costly dynamics memory allocations. Experiments with a fluid simulator code show that the gain for sorting can be significant, even if the overall performance benefit is limited due to the reduced weight of sorting.

The mechanism could benefit as well to other particle-based simulations such as crowds, molecular dynamics, granular materials, or boids simulations. Future work also include evaluating the benefits of PMAs for other particle based simulations, and developing parallel PMA implementations for multi-core CPUs and GPUs.

Figure 1 :

 1 Figure 1: Particles (in brown) placed into a collision detection grid in 2D. Cells are Z-curve ordered (blue segmented line).

  [START_REF] Bender | Cache-oblivious b-trees[END_REF][START_REF] Bender | An adaptive packedmemory array[END_REF]. A PMA consists in an array of size N . Gaps are empty array elements, opposed to filled array elements, called elements, that actually contain the sorted data. The PMA is evenly divided into c The Eurographics Association 2012.

Figure 2 :

 2 Figure 2: PMA representation with S = 4 segments of size 5. Every level of the segments tree owns its minimum and maximum element densities.

Figure 4 :

 4 Figure4: Collision detection grid in 2D with particles (in brown) sorted according to their cell index (blue Zcurve) in two dense arrays (top right) or using two PMA data structures,one for storing each non empty cell and a pointer toward the cell first particle, and one for the particles (bottom right).

c

  The Eurographics Association 2012.

  Figure6: Total execution times (s) for a CBD simulation of 180 000 particles with native code when particles are never sorted in memory after initialization (No Sort), when they are sorted every 100 time step (Sort Every 100), every time step (Reference), and using the PMA. On the left, we show the execution time for 1 000 time steps. On the right, we focus on the 580 -720 range to highlight differences between sort politics.

  element in the 2 l segments window This scheme ensures that any insertion or deletion is performed in O(log 2 N ) amortized element moves and that for special patterns, like random insertions, the cost is reduced to O(log N )[START_REF] Bender | Cache-oblivious b-trees[END_REF][START_REF] Bender | An adaptive packedmemory array[END_REF].

	4)
	Algorithm 2 Pseudo algorithm for one element in-
	sertion inside a PMA.
	binary search the segment s the new element should
	go into
	if τ0 exceeded then
	search for the l window that fits τ l
	if τ h not respected then
	double PMA, l ← h
	rebalance elements
	else
	rebalance window of 2 l segments
	end if
	end if
	search and insert

  Table 1 gives average execution times for 5% moves on various number of elements. The PMA outperforms the quicksort

c The Eurographics Association 2012.

Table 3 :

 3 Impact of the cell size (size), and of the number of neighbor cells tested (#nb) on the number of tested pairs and the collision detection (CD) time (ms) with native Fluids implementation.

	size, #nb tested pairs valid pairs ratio CD
	2h, 8	36878962	2691708	13.7	298
	h, 27	16472058	2691708	6.1	232
	h, 14	8236029	2691708	3.1	145

Table 4 :

 4 Relative weights of the different steps of the simulation. % values are averaged over 4000 time steps simulating 180 000 particles with Fluids code sorting particles at every time step.

		1
						Reference
	Execution time [s]	0.2 0.4 0.6 0.8				PMA % moves	30% 10% 20%	% moves
		0	0	100	200	300 0%
				Frames
	Figure 5: Sort execution time (s) and move ratio (%)
	per iteration during a CBD simulation of 2 900 000
	particles, for Fluids with particles sort at every time
	step (Reference) and for the PMA.
			step			%
			sort			4.5
			detect neighbors	49.1
			compute pressures	6.3
			compute forces	36.8
			update positions		3.2

c The Eurographics Association 2012.

Acknowledgments

We would like to warmly thank Swann Perarnau for his contribution to code development. This work was partly funded by ANR projet Repdyn 09-COSI-011-05.