Dominik M Wittmann

Carsten Marr

Fabian J Theis

Biologically meaningful update rules increase the critical connectivity of generalized Kauffman networks

Keywords: dynamic systems on graphs, multilevel logic, thermodynamics, transcriptional gene regulation 2010 MSC: 80A30, 82C20, 82C27, 92C42

We generalize random Boolean networks by softening the hard binary discretization into multiple discrete states. These multistate networks are generic models of gene regulatory networks, where each gene is known to assume a finite number of functionally different expression levels. We analytically determine the critical connectivity that separates the biologically unfavorable frozen and chaotic regimes. This connectivity is inversely proportional to a parameter which measures the heterogeneity of the update rules. Interestingly, the latter does not necessarily increase with the mean number of discrete states per node. Still, allowing for multiple states decreases the critical connectivity as compared to random Boolean networks, and thus leads to biologically unrealistic situations.

Therefore, we study two approaches to increase the critical connectivity.

First, we demonstrate that each network can be kept in its frozen regime by sufficiently biasing the update rules. Second, we restrict the randomly chosen update rules to a subclass of biologically more meaningful functions. These functions are characterized based on a thermodynamic model of gene regu-

Introduction

Boolean networks are a class of discrete dynamical systems. The system's N variables take only discrete values 0 or 1 and develop in discrete time steps.

At each time point, the value of a variable is determined by a so-called update rule that deterministically depends upon the values of some of the other variables, the so-called inputs, at the previous time point. [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF] proposed random Boolean networks -so called Kauffman networks (KN)as generic models for large-scale gene regulatory networks [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF].

Computational experiments showed that these networks exhibit surprisingly ordered structures and are able to give insights into biological phenomena such as cell replication or lineage differentiation.

Interest in KNs was rekindled as their close relation to classical models from statistical mechanics was realized. In a number of studies [START_REF] Derrida | Random networks of automata: a simple annealed approximation[END_REF][START_REF] Derrida | Phase transitions in two dimensional Kauffman cellular automata[END_REF][START_REF] Flyvbjerg | An order parameter for networks of automata[END_REF] the selforganizing capacity of KNs was analyzed. It was shown that depending on their connectivity, KNs exhibit ordered (frozen) as well as chaotic behaviors with a critical boundary separating both regimes. The ordered regime is characterized by small stable attractors, whereas in the chaotic regime long-periodic orbits frequently occur. These properties render both regimes unfavorable for the evolution of living organisms. Consequently, Kauffman promoted the idea of "living at the edge of chaos" [START_REF] Kauffman | The origins of order: Self organization and selection in evolution[END_REF]. Interestingly, the critical connectivity of KNs is 2, which agrees well with the average connectivities of well-studied gene regulatory networks, e.g. in E. coli, S. cerevisiae and B. subtilis [START_REF] Balleza | Critical dynamics in genetic regulatory networks: examples from four kingdoms[END_REF].

Besides studies of large-scale random Boolean networks, many smalland medium-scale Boolean models of different biological processes have been manually curated and analyzed [START_REF] Albert | The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in drosophila melanogaster[END_REF][START_REF] Li | The yeast cell-cycle network is robustly designed[END_REF][START_REF] Saez-Rodriguez | A logical model provides insights into t cell receptor signaling[END_REF]Wittmann et al., 2009a). Here however, the discretization of continuous biological quantities, such as mRNA or protein concentrations, into binary 'on'-'off' categories is often arbitrary as well as insufficient. In fact, it has been demonstrated that genes may well have more than two functionally different expression levels [START_REF] Setty | Detailed map of a cis-regulatory input function[END_REF].

To alleviate this issue, in some modeling applications the Boolean categories are extended to multiple discrete states, leading to multistate logical models [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF][START_REF] Sánchez | A logical analysis of the Drosophila gap-gene system[END_REF][START_REF] Fauré | Modular logical modelling of the budding yeast cell cycle[END_REF].

Studies of KNs with multiple states, so-called multistate Kauffman networks (MKN), are scarce [START_REF] Solé | Phase transitions in random networks with multiple states[END_REF] and restricted to the biologically implausible case where all nodes have the same number of discrete states.

In this case it has been observed that the critical connectivity tends to 1, as the number of states grows. Hence, it drops below the average connectivities we typically find in gene regulatory networks, cf. again [START_REF] Balleza | Critical dynamics in genetic regulatory networks: examples from four kingdoms[END_REF], leading to biologically unrealistic scenarios. For (Boolean) KNs two approaches are known that can lead to an increased critical connectivity.

(i) First, it has been shown that for any connectivity a KN can be kept in the frozen regime by biasing the update rules towards one of the discrete states [START_REF] Bastolla | Closing probabilities in the Kauffman model: an annealed computation[END_REF].

(ii) Second, the issue of biologically meaningful update rules has attracted considerable attention, for a brief review see section 6.3.

In this contribution, we generalize these approaches to MKNs and ask how they affect their critical connectivity. We begin by investigating a very general class of MKNs in section 2. In particular, contrary to previous studies, the number of discrete states is not fixed for all nodes but follows a distribution. We show that the above mentioned transition between a frozen and a chaotic regime takes place in MKNs of this class. The critical boundary is analytically determined as a relation between the mean connectivity and a parameter describing the heterogeneity of the update rules. Our results are an extension of the (binary) Boolean case. They show that allowing nodes to assume multiple states lowers the critical connectivity as compared to (Boolean) KNs, the limit being 1. Interestingly however, the critical connectivity does not necessarily decrease with the mean number of states per node. With respect to our motivating question, we find in section 3 thatsimilar to (Boolean) KNs -each MKN can be kept in the frozen regime by biasing the update rules towards one of the discrete states: Approach (i) can be generalized to MKNs.

The main part of this manuscript is devoted to the generalization of approach (ii). In section 4 we propose a tentative characterization of biologically plausible update functions for MKNs. This characterization is motivated by a common property of the update rules from some manually curated multistate logical models, which we call single-switch condition. It is further corroborated by a thermodynamic model of transcriptional gene regulation.

Using this model we show that under certain assumptions regarding cooperative effects the promotor activity as a function of the transcription factor concentrations always satisfies the single-switch condition. From a general point of view, our thermodynamic considerations link discrete and continuous models of gene regulatory networks.

The effect of restricting update rules to the class of biologically more meaningful single-switch functions is studied in section 5. Analytic results demonstrate that the use of single-switch functions indeed increases the critical connectivity in MKNs as compared to fully random update rules. Hence, our motivating question can be answered positively also with respect to approach (ii). The usage of single-switch functions causes further interesting effects. We present an example, in which for an increasing number of states the critical connectivity increases if single-switch functions are used, while it decreases if the update rules are chosen randomly. Section 6 concludes our manuscript by a non-technical summary of our results and their discussion, in particular, in the light of previous work.

Phase transitions in multistate Kauffman networks

In this section we introduce a very general class of MKNs and investigate it with respect to critical phenomena.

A general class of multistate Kauffman networks

A discrete multistate network consists of N variables x 1 , x 2 , . . . , x N each taking values in a discrete finite set Σ i with cardinality S i , w.l.o.g. Σ i = {0, 1, . . . , S i -1}, i = 1, 2, . . . , N. Time is discretized, t = 0, 1, . . ., and we denote the variables' time courses by x i (t). The i-th variable is influenced by K i inputs x i1 , x i2 , . . . , x iK i and its value at time t + 1 is given by an update

function f i : K i k=1 Σ ik → Σ i , x i (t + 1) = f i (x i1 (t), x i2 (t), . . . , x iK i (t)) .
We call K i the connectivity of x i . The state of a multistate network at time t is denoted by

X(t) = (x 1 (t), x 2 (t), . . . , x N (t)) ∈ N i=1 Σ i .
A Boolean network is the special case of a multistate network where all S i = 2.

We now define MKNs generalizing the definition of KNs as given e.g. in [START_REF] Aldana | Boolean dynamics with random couplings[END_REF]. A MKN is a multistate network where (K1) the K i are chosen randomly from a probability distribution P in (K), K = 1, 2, . . . , K max , K max ≤ N , (K2) the K i inputs x i1 , x i2 , . . . , x iK i of x i are chosen randomly with uniform probability from among the system's variables x 1 , x 2 , . . . , x N , (K3) the number of states S i of x i is chosen randomly from a probability distribution P nos (S), S = 2, 3, . . . , S max , (K4) the values of f i are independent and chosen randomly from a probability distribution P S i (s), s ∈ Σ i .

Note that in (K4) the distribution P S i does not depend on the node x i but only on the number of states S i of node x i . In particular, in the Boolean case the update rules evaluate to 0 with a certain probability w and to 1 with probability 1w.

The critical boundary

From (K4) it follows that the probability p S i for the function f i to yield two different values for two different arguments depends only on P S i and is given by

p S i = s∈Σ i P S i (s) (1 -P S i (s)) .
(1)

We remark that p S i becomes maximal if P S i is the discrete uniform distribution,

P S i (s) = 1 S i , s ∈ Σ i , (2)
cf. Appendix A. In this case,

p S i = S i -1 S i ,
which also goes to show that p S i < 1 for any choice of P S i .

Following e.g. [START_REF] Aldana | Boolean dynamics of networks with scale-free topology[END_REF] we distinguish between the frozen and the chaotic regime by looking at the overlap o(t) between two distinct time courses

X(t) = (x 1 (t), x 2 (t), . . . , x N (t)) and X (t) = (x 1 (t), x 2 (t), . . . , x N (t))
of the variables x 1 , x 2 , . . . , x N . The overlap at time t is defined as the fraction of variables that take the same value in X(t) and X (t), i.e.

o(t) def = 1 N N i=1 δ x i (t),x i (t) . ( 3
)
We now express o(1) in terms of o(0). Fixing some i, there are two (mutually exclusive) possibilities that at t = 1 we have x i (1) = x i (1):

1. All K i inputs are equal at t = 0. This occurs with probability o(0) K i .

2. Not all inputs are equal at t = 0 but f i still yields the same value. This occurs with probability 1o(0

) K i (1 -p S i ).
Hence, the probability for "x i (1) = x i (1)" is given by o(0

) K i + 1 -o(0) K i (1 -p S i )
and depends only on K i and p S i . It follows that

o(1) = Kmax K=1 P in (K) Smax S=2 P nos (S) o(0) K + 1 -o(0) K (1 -p S ) .
We now apply a trick devised by [START_REF] Derrida | Random networks of automata: a simple annealed approximation[END_REF] called annealed approximation. Although the model is quenched, i.e. fixed for all time points, we assume that after each time step a new model is generated according to (K1)-(K4). It has been shown that in the thermodynamic limit N → ∞ the quenched and the annealed model show the same phase transition [START_REF] Hilhorst | On the approach of the stationary state in Kauffman's random Boolean network[END_REF]. Therefore, for N K max we can write

o(t + 1) = O (o(t)) def = Kmax K=1 P in (K) Smax S=2 P nos (S) o(t) K + 1 -o(t) K (1 -p S ) (4)
at each time step.

In the limit t → ∞ this dynamic equation becomes a fixed point equation

o = O(o)
for the stationary value of the overlap. Note that o = 1 is always a fixed point. If o is attractive, initial differences between time courses vanish over time and the system is in an ordered regime, otherwise it is chaotic.

Hence, the ordered and chaotic regimes as well as the critical boundary can be characterized by a stability analysis of o . As detailed in Appendix B this yields

pK ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 1 ordered regime = 1 critical boundary > 1 chaotic regime , (5)
where

p def = Smax S=2 P nos (S) p S ( 6
)
and

K def = Kmax K=1 P in (K) K .
In other words, the critical connectivity is given by

K crit = 1 p . ( 7
)
We observe that the phase transition depends only on the first moments K and p. Typically, K can be easily computed once the underlying network structure is known. The computation of p is more involved.

The situation of fully random update rules without bias

As a preparation for section 3 let us now study the multistate counterpart of the original KN as introduced by [START_REF] Kauffman | Metabolic stability and epigenesis in randomly constructed genetic nets[END_REF]. This is to say we consider the completely unbiased situation where the P S i are the discrete uniform distributions from (2). We already remarked that in this case

p S i = S i -1 S i . ( 8
)
Substituting ( 8) in (6) yields

p = Smax S=2 P nos (S) S -1 S = 1 -μ ,
where

μ def = Smax S=2 P nos (S) 1 S .
Clearly, for any P nos , μ is bounded, 0 < μ ≤ 1/2. For the critical connectivity from ( 7)

K crit = 1 1 -μ , ( 9
)
this implies 1 < K crit ≤ 2. Moreover, K crit lowers down to 1 as μ becomes small.

Hence, the critical connectivity in unbiased MKNs is smaller than in unbiased (Boolean) KNs. This important result serves as the motivation for our investigations in sections 3, 4 and 5. Before, let us visualize this result in two examples.

Examples

Let us conclude this section by computing μ in two special cases of P nos .

• First, a delta distribution

P nos (S) = δ S,S , S = 2, 3, . . . , S max , ( 10
)
for a fixed number of states S per node. In this case, μ = 1/S, and from (9) we obtain the condition

K crit = S S -1 , ( 11
)
which was already found in [START_REF] Solé | Phase transitions in random networks with multiple states[END_REF]. Note that in the Boolean case, S = 2, we recover the critical connectivity 2 of KNs. For an increasing (fixed) number of states S the critical connectivity approaches 1.

• Second, a discrete uniform distribution

P nos (S) = 1 S max -1
, S = 2, 3, . . . , S max .

In this case,

μ = Smax S=2 1 S max -1 1 S = H Smax -1 S max -1 ,
where H Smax is the S max -th harmonic number, and from ( 9) we obtain the criticality condition

K crit = S max -1 S max -H Smax . ( 12
)
It holds that K crit → 1 as S max → ∞. A better idea of this convergence is given by the well-known approximation

H Smax = ln (S max ) + γ + O(S -1 max ) ,
where γ is the Euler-Mascheroni constant γ ≈ 0.5772.

The critical connectivities from (11) and ( 12) are plotted in Figure 1A.

We observe that in both examples a growing (fixed) number of states S or a growing upper bound S max , respectively, lower the critical connectivity down to 1. In the following sections we generalize the approaches (i) and (ii)

from section 1 to MKNs and investigate if they lead to an increased critical connectivity.

Multistate Kauffman networks with biased update rules

In the original KN, each Booelan update function was drawn randomly from the ensemble of all Boolean functions according to the uniform distribution. In subsequent studies, however, it is mostly assumed that each function evaluates to 0 with a certain probability 0 < w < 1 (and to 1 with probability 1w). Biologically speaking, it is assumed that a gene is expressed with probability w. A natural generalization hereof to multistate systems is the distribution

w c ri t 1 ( K , 2 ) w c ri t 2 ( K , 2 ) w c ri t 1 ( K , 3 ) w c ri t 2 ( K , 3 ) w c ri t ( K , 3 )
P S i (s) = ⎧ ⎨ ⎩ w s= 0 1 -w S i -1 s > 0 (13)
for some 0 < w < 1.

In the Boolean case, it turned out that for each connectivity a system can be kept in the frozen regime by choosing w sufficiently large or small [START_REF] Bastolla | Closing probabilities in the Kauffman model: an annealed computation[END_REF]. Is this still true for MKNs?

The critical boundary

To answer this question, let us assume that each P S i is given by ( 13) for some 0 < w < 1. Then (1) becomes

p w S i = 2w(1 -w) + (1 -w) 2 S i -2 S i -1 ,
and substituting this in (6) gives

p w = Smax S=2 P nos (S) 2w(1 -w) + (1 -w) 2 S -2 S -1 = 2w(1 -w) + (1 -w) 2 (1 -μ -1 ) , ( 14
)
where

μ -1 def = Smax S=2 P nos (S) 1 S -1 .
We observe that 0 < μ -1 ≤ 1, where the upper bound is exact iff P nos (S) = δ S,2 , i.e. in the Boolean case.

Substituting ( 14) in ( 5) gives rise to a quadratic equation in w. Solving this equation for w crit (K, μ -1 ) yields (cf. Appendix C):

(a) For K < 1 + μ -1 the system is always frozen.

(b) For K = 1 + μ -1 the system is critical if and only if

w = w crit (1 + μ -1 , μ -1 ) = μ -1 μ -1 + 1 ,
otherwise it is frozen.

(c) In the Boolean case, i.e. μ -1 = 1, and for K > 1 + μ -1 = 2 the critical boundary is described by

w crit 1 2 K, 1 = 1 2 1 ∓ 1 - 2 K ∈ (0, 1) .
For

w crit 1 K, 1 < w < w crit 2 K, 1
the system is chaotic, otherwise it is frozen. This agrees with previous results about Boolean KNs, see e.g. [START_REF] Bastolla | Closing probabilities in the Kauffman model: an annealed computation[END_REF].

(d) For μ -1 < 1 and 1/(1μ -1 ) > K > 1 + μ -1 the critical boundary is described by

w crit 1 2 K, μ -1 = 1 μ -1 + 1 μ -1 ∓ 1 - μ -1 + 1 K ∈ (0, 1) .
For

w crit 1 K, μ -1 < w < w crit 2 K, μ -1 the system is chaotic, otherwise it is frozen. For K = 1/(1 -μ -1 ),
w crit 1 = 0 is not a valid solution for w anymore, as we require w ∈ (0, 1).

(e) For μ -1 < 1 and K ≥ 1/(1μ -1 ) the critical boundary is described by

w crit K, μ -1 = 1 μ -1 + 1 μ -1 + 1 - μ -1 + 1 K ∈ (0, 1) ,
for larger values of w the system is frozen, otherwise it is chaotic. It can easily be seen that w crit (K, μ -1 ) is monotonous in both its arguments and that w crit → 1 as K → ∞.

We chose to solve the criticality condition for w (and not K crit ), as from this presentation we easily see that for any values of K and μ -1 a system is kept in its frozen regime by choosing a sufficiently large w. Thus, our motivating question can be answered in the affirmative: Biasing update rules towards one of the discrete states indeed increases the critical connectivity, from a theoretical point of view, even beyond any bound.

Examples and visualization

For clarification and illustration purposes, let us again consider the example where P nos follows the delta distribution from (10) for a fixed number of states S. This entails μ -1 = 1/(S -1), and we can express the classification from section 3.1 in terms of S, cf. Appendix D. This classification is visualized in Figure 1B. We observe that w crit (K, S) is monotonously increasing in S and K, and that w crit → 1 as K → ∞.

To further visualize our results, we generated a frozen, critical and chaotic MKN, for details see caption of

A class of biologically meaningful update rules

In the last sections, the update functions f i were chosen from the set of all possible update rules according to (K4). However, as argued below, especially for higher number of states most of these update rules are arguably not plausible from a biological point of view. In this section, we characterize a class of biologically more meaningful update rules for MKNs. We illustrate the main characteristics of this class in section 4.1. Subsequently, in section 4.2, we present a more formal derivation using a thermodynamic model of gene regulation. MKNs with biologically meaningful update rules are then investigated in section 5, especially with respect to their critical connectivity. this question let us take a look at the manually curated update rules from multistate logical models, e.g. of the yeast cell cycle [START_REF] Fauré | Modular logical modelling of the budding yeast cell cycle[END_REF] or the Drosophila Gap-gene System [START_REF] Sánchez | A logical analysis of the Drosophila gap-gene system[END_REF]. The "least common denominator" of these rules is, arguably, that there is only one threshold for each input at which the function changes its output. In other words, there is only a single switching point along each input dimension.

Update rules with single-switch inputs

We will call functions satisfying this condition single-switch functions. More formally, a single-switch function f with inputs

x 1 ∈ Σ 1 , x 2 ∈ Σ 2 , . . . , x K ∈ Σ K has the property that for each 1 ≤ k ≤ K there exists a threshold θ k ∈ Σ k such that sgn (x k -θ k ) = sgn (x k -θ k ) implies f (x 1 , . . . , x k , . . . , x K ) = f (x 1 , . . . , x k , . . . , x K ).
An example of a single-switch function is shown in Figure 2B.

One might wonder why we allow multiple states for the input variables if in the end they behave essentially like Boolean (binary) variables. This is certainly a valid argument if we consider only one update rule. Typically however, a transcription factor influences more than one gene and the thresholds for these regulations can differ. This is, in fact, precisely the reasoning that led to the development of multistate logical models. We would like to refer the reader to the review by [START_REF] Thomas | Regulatory networks seen as asynchronous automata: A logical description[END_REF] where the author addresses inter alia the question "When and How Should One Use Variables with More Than Two Values?"

We argue that imposing the single-switch condition is a reasonable restriction of update rules and that MKNs with single-switch update rules are closer to biological reality than MKNs with fully random update rules, see section 6.3 for further discussion. Therefore, we believe it worthwhile to investigate MKNs with single-switch update rules and we will do so in section 5. Before, let us further and more formally validate the single-switch criterion.

A thermodynamic model of gene regulation

In order to provide a solid theoretical foundation for our characterization of meaningful update rules, we now describe a thermodynamic model of gene regulation, leaning heavily on [START_REF] Sneppen | Physics in molecular biology[END_REF]. Here, we present only the most important results and formulas. More detailed computations and derivations can be found in Appendix E.

Let us consider a prototypic gene x that is regulated by K regulators

x 1 , x 2 , . . . , x K . For each regulator x k , 1 ≤ k ≤ K, there are n k binding sites in the promotor region of x, cf. Figure 2G. We assume maximal homotypic cooperativity, which is to say, we assume that each transcription factor x k , 1 ≤ k ≤ K, can bind only as a polymer consisting of n k monomers. Moreover, we neglect heterotypic cooperativity, i.e. cooperative effects between two different transcription factors. For further discussion of these assumptions see section 6.4. We thus have 2 K accessible promotor configurations, which we denote by vectors ω ∈ {0, 1} K , the k-th component ω(k) indicating whether or not the k-th regulator polymer is bound.

As detailed in Appendix E we can express the probability of a certain promotor configuration ω as

π ω = k|ω(k)=1 h + ([x k ] , n k , A k ) • k|ω(k)=0 h -([x k ] , n k , A k ) , ( 15
)
where

h + ([x k ] , n k , A k ) = A k [x k ] n k 1 + A k [x k ] n k and h -([x k ] , n k , A k ) = 1 1 + A k [x k ] n k
ω, r ω ) ∈ R K+1 , B I (rω) (y 1 , y 2 , . . . , y K ) = ω r ω k|ω(k)=1 y k • k|ω(k)=0 (1 -y k ) . ( 16
)
Using ( 15) and recalling that h

-([x k ] , n k , A k ) = 1 -h + ([x k ] , n k , A k ), we can
express the overall transcription rate r as

r = ω π ω r ω = B I (rω) (h + ([x 1 ] , n 1 , A 1 ) , h + ([x 2 ] , n 2 , A 2 ) , . . . , h + ([x K ] , n K , A K )) , (17)
i.e. as the composition of the multilinear interpolation with Hill functions.

Hence, we can nicely separate linearities and non-linearities. Non-linearities appear only as single-component non-linearities and take the form of standard Hill functions. The different inputs are multilinearly coupled.

In discrete models of gene regulation the overall transcription rate r of gene x as a (continuous) function of the transcription factor concentrations

[x k ], 1 ≤ k ≤ K,
is mimicked by a discrete update rule. This process is formalized in the next section. It will turn out that update rules obtained from transcription rates r as in ( 17) by such a discretization are always single-switch functions.

The transition from continuous to discrete models of gene regulation

As a first step towards a discrete model we replace the sigmoidal Hill 17) by Heaviside step functions

functions h + ([x k ] , n k , A k ) in (
h step [x k ] , D micro k = ⎧ ⎨ ⎩ 0, [x k ] < D micro k 1, [x k ] ≥ D micro k , ( 18
)
where

D micro k = n k 1/A k is the so-called microscopic dissociation constant,
which indicates the concentration level for half-maximal activation, cf. Figure 2C. This results in piecewise linear functions

B P W (rω) = B I (rω) h step [x 1 ] , D micro 1 , . . . , h step [x K ] , D micro K ( 19
)
as introduced by Glass and [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF].

We can now easily switch to discrete variables x ∈ Σ = {0, 1, . . . , S -1} and x k ∈ Σ k = {0, 1, . . . , S k -1} by choosing transcription rates r ω ∈ Σ and defining our update function f of x via f :

K k=1 Σ k -→ Σ (x 1 , x 2 , . . . , x K ) -→ B P W (rω) (x 1 , x 2 , . . . , x K ) . ( 20
)
In the discrete case we can w.l.o.g. choose the thresholds

D micro k from Σ k .
Further note that in definition (18) the threshold value itself is contained in the preimage of 1. Hence, in the discrete case, the set of non-degenerate thresholds is Σ k \ {0}, as for these values there is a non-empty left-as well as right-hand side.

The resulting f from ( 20) is clearly a single-switch function (with switching point D micro k along the x k -dimension) as defined above in section 4.1.

Thus, we could corroborate our characterization of meaningful update rules.

In the following, we will consider update rules satisfying the single-switch condition. In particular, we ask if their usage leads to an increased critical connectivity in MKNs. Before, let us further illustrate our derivation of single-switch functions using a simple example.

Example

Let us show how the single-switch function shown in Figure 2B could have arisen from the thermodynamic model described above. Assume that a gene x is regulated by two transcription factors x 1 and x 2 , which can bind to the promotor region of x as n 1 = 3-mers and n 2 = 4-mers, respectively, (maximal homotypic cooperativity). The bonding affinity between x 1 (x 2 )

and the promotor region determines the threshold

D micro 1 (D micro 2
) for the regulation of x by x 1 (x 2 ). Here, we take D micro 1 = 1/3 and D micro 2 = 1/2 in unit-normalized concentrations. We denote the possible promotor configurations by (0, 0), (1, 0), (0, 1) and (1, 1), the first (second) component indicating whether or not x 1 (x 2 ) is bound. We choose (unit normalized) transcription rates r (0,0) = 0, r (1,0) = 0.25, r (0,1) = 0.75, r (1,1) = 1 for these different configurations.

Figure 2D shows the transcription rates r ω , ω ∈ {0, 1} 2 , and their mul- would lie between the first and the second. Finally, assume that gene x is also a transcription factor and regulates four genes at different thresholds. This leads to 5 discrete states for x. After this discretization, the piecewise linear function from Figure 2F becomes the single-switch function shown in Figure 2B.

Multistate Kauffman networks with single-switch functions

In the following, we no longer draw the update rules randomly from the set of all possible rules, but modify (K4) to (K4') Each f i is a single-switch function with non-degenerate thresholds randomly chosen from Σ ik \ {0} and output values (transcription rates) chosen according to the distribution P S i .

Does this increase the critical connectivity as compared to ( 7)? To answer this question let us analytically describe the critical boundary.

The critical boundary

Due to the restricted architecture of the update rules, we have to slightly modify our reasoning from section 2.2. We begin by deriving an appropriate iteration for o(t) from (3), which was defined as the overlap between two time courses X(t) and X (t). To this end, let us consider the i-th node, 1 ≤ i ≤ N , and its k-th input, 1 ≤ k ≤ K i . We denote the probability that for two different values x ik = x ik of this input we have

h step (x ik , θ) = h step (x ik , θ),
where θ ∈ Σ ik \ {0} is a randomly chosen threshold, by q S ik . For a way to compute q S ik , see Appendix F. Let

q def = Smax S=2 P nos (S) q S .
Clearly, q = 1 iff for all S with P nos (S) = 0 the support of P S consists of exactly two points, i.e. iff the network is essentially a Boolean network.

Analogously to section 2.2, we reason that, given o(t), there are two (mutually exclusive) possibilities that at t + 1 we have x i (t + 1) = x i (t + 1).

1. For each of the K i inputs either of the following is true:

1.1. The input is contained in the overlap at t, i.e. x ik (t) = x ik (t). This occurs with probability o(t).

1.2. The input is not contained in the overlap, x ik (t) = x ik (t), but still h step (x ik (t), θ) = h step (x ik (t), θ). This occurs with probability

(1 -o(t)) (1 -q S ik ).
The overall probability for possibility 1 is given by

Smax S=2 P nos (S) [o(t) + (1 -o(t)) (1 -q S )] K i = [o(t) + (1 -o(t)) (1 -q)] K i .
2. Possibility 1 is wrong, but f i still yields the same value. Using (1), this occurs with probability 1

-[o(t) + (1 -o(t)) (1 -q)] K i (1 -p S i ).
We can now define an iteration for the overlap o(t).

o(t + 1) = O switch (o(t)) def = Kmax K=1 P in (K) Smax S=2 P nos (S) [o(t) + (1 -o(t)) (1 -q)] K + 1 -[o(t) + (1 -o(t)) (1 -q)] K (1 -p S ) .
Clearly, o = 1 is a fixed point. A stability analysis of this fixed point very similar to Appendix B yields the following phase transition:

pKq ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 1 ordered regime = 1 critical boundary > 1 chaotic regime .
Consequently, the critical connectivity is given by

K crit = 1 p q . ( 21
)
Comparing ( 21) to (7) we observe that the restriction of update rules to single-switch functions introduces the factor 1/q into the criticality condition.

Since q < 1 in the non-Boolean case, this leads to an increased critical connectivity for MKNs, which answers our initiatory question.

Example

Let us conclude by considering the special case in which P nos is the delta distribution from (10) for a fixed number of states S, and P S is the discrete uniform distribution from (2). The (mean) probability q = q S then computes to q = 1 3

S + 1 S -1 , cf. Appendix F.
Dividing the critical connectivity from (11) by q yields

K crit = 3 S S + 1 . ( 22
)
In Figure 2H the critical connectivity from ( 22) and -for comparisonthe critical connectivity from ( 11) are shown. For larger values of S the use of single-switch functions leads to a threefold increase in the critical connectivity as compared to fully random update rules. Interestingly, while for fully random update rules the critical connectivity decreases as S → ∞, it increases if single-switch functions are used.

Summary and discussion

KNs as generic models of gene regulatory networks are well studied and able to explain certain aspects of biological processes such as cell replication or lineage differentiation. Still, in many respects these models are arguably too crude a simplification. In this contribution, we proposed modifications of the original KN in order to enhance its biological plausibility.

Summary

MKNs can be seen as a natural generalization of (Boolean) KNs and allow us to soften the hard binary discretization. We investigated a very general class of MKNs and demonstrated that they exhibit a phase transition from frozen to chaotic behavior. In the frozen regime, networks are robust in the sense that they tolerate perturbations. In the chaotic regime, even small perturbations will lead to the divergence of trajectories and thus to qualitatively different behaviors. The critical boundary between both regimes was determined analytically. In its most general representation (7) the critical connectivity is inversely proportional to the mean heterogeneity parameter p from (6) of the update rules. This parameter becomes maximal if the update rules are unbiased, i.e. assume each value in the state space with equal probability. It is minimal (zero), if the update rules are constant functions.

A nice intuitive interpretation of the phase transition ( 5) is given in terms of "damage spreading". Assume that at time t the state of a MKN is damaged in one node, i.e. the state of this node is altered. On average this node affects K other nodes at the next time step t + 1. The parameter p is the mean probability that a change of input indeed leads to a change of output of an update rule. The product Kp thus gives the mean number of damaged (changed) nodes at time t + 1. In the case Kp > 1 damage will spread through the network, in the case Kp < 1 the network is able to "repair" damage over time.

Our results show that MKNs have smaller critical connectivities than the original (Boolean) KNs. Yet we also demonstrated that each system can be kept in its frozen regime by putting a sufficiently heavy bias on one of the states, cf. Figure 1B. In the Boolean case, the critical weights w crit 1/2 are, of course, symmetric about 1/2. In MKNs there still is a range of connectivities where we can freeze a network by choosing either sufficiently large or sufficiently small weights. For higher connectivities, however, the option of choosing small weights ceases to exist. Intuitively speaking, even if we set the weight to zero, the heterogeneity among the remaining states would still be too large. Mathematically, this is reflected in the solution w crit 1 becoming negative. From a biological point of view, this may indicate that in real genetic networks the update functions have a base level of activation and deviations thereof constitute well-defined exceptions.

In a next step, we restricted the randomly chosen update rules to a subclass of biologically more meaningful functions, so-called single-switch functions. These were characterized based on a thermodynamic model of gene regulation. The crucial property of single-switch functions is that their output value changes only once along each input dimension. Hence, there is a chance that a single-switch function does not "realize" a change in one of its inputs, simply because the input remains on the same side of its switching point. In terms of "damage spreading" this implies that damage is less likely to spread. (We remark that this is precisely what is measured by the additional factor 1/q in (21).) Consequently, we expect MKNs with single-switch update rules to be more robust. Indeed we found that the use of single-switch functions increases the critical connectivity; in a specific example by a factor of up to three. Similar to the study of nested canalyzing functions [START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF], our results show the importance of the update functions' architecture for the dynamic behavior of KNs.

Previous work on multistate Kauffman networks

When employing multistate logical models one needs to make a fundamental decision about how update rules act on variables. Either the value of an update rule determines the absolute value of a variable at the next time step or it determines the change of the variable relative to the previous state, i.e. whether the variable is assigned the next higher or the next lower state.

KNs with multiple states using the latter update policy are called random walk networks. Analyses of these networks revealed a phase transition be-tween a chaotic and an ordered regime similar to the standard KN [START_REF] Luque | Random walk networks[END_REF][START_REF] Ballesteros | Order-disorder phase transition in randomwalk networks[END_REF]. In the chaotic regime, variables follow random-walk like trajectories, which gives the name to these networks. In here, we did not consider random walk networks but KNs with multiple states using the first update policy; we referred to them as MKNs.

Previous studies [START_REF] Solé | Phase transitions in random networks with multiple states[END_REF] of MKNs are restricted to the biologically unreasonable case in which all nodes have the same number of states.

In agreement with our more general results, it has been shown that in this case an increasing number of states decreases the critical connectivity down to 1. It is tempting to assume that this stays true (at least in a qualitative sense) when replacing the (fixed) number of states by the mean of a nondegenerate distribution. However, our results demonstrate that, in general, this is wrong. From the critical connectivity in ( 9), for instance, we deduce that here the crucial parameter is not the mean number of states, but the mean of the reciprocal number of states μ.

Generally speaking, the mean of a strictly positive, non-constant random variable R and the mean of 1/R are related by

E 1 R > 1 E (R) ,
for a proof see e.g. [START_REF] Kendall | Kendall's advanced theory of statistics[END_REF]. Furthermore, we can easily come up with a distribution

P (R = 2) = r -1 r and P R = r 2 = 1/r for r ∈ N, such that E(R) → ∞ but E (1/R) → 1/2 as r → ∞.
Setting P nos = P , the critical connectivity from (9) increases up to 2 as r → ∞, although the mean number of states grows beyond any bound.

Biologically meaningful update functions

The issue of biologically meaningful update rules for Boolean networks has attracted considerable attention. In [START_REF] Harris | A model of transcriptional regulatory networks based on biases in the observed regulation rules[END_REF] it has been observed that in lower organisms most genes have one so-called canalyzing (forcing) input, such that for one value of this input, the output value is fixed.

Only if this input is not canalyzing, do the other inputs become relevant. In [START_REF] Kauffman | Random Boolean network models and the yeast transcriptional network[END_REF] the concept of canalyzing functions was extended to nested canalyzing functions. This was done mostly for technical reasons as for nested canalyzing rules simple sampling schemes exist. In a subsequent study, networks with nested canalyzing rules were shown to be always stable [START_REF] Kauffman | Genetic networks with canalyzing Boolean rules are always stable[END_REF].

A different class of biologically plausible Boolean rules was defined in [START_REF] Raeymaekers | Dynamics of Boolean Networks Controlled by Biologically Meaningful Functions[END_REF] based on the assumption that each transcription factor is either an overall activator or an overall inhibitor. If one furthermore assumes that the activatory and inhibitory regulators of a gene interact in an additive fashion, one can model genetic networks by so-called random threshold networks [START_REF] Rohlf | Criticality in random threshold networks: annealed approximation and beyond[END_REF][START_REF] Szejka | The phase diagram of random threshold networks[END_REF]. Here, the update function of a node is a weighted sum over its inputs (with positive and negative weights depending on whether the input is an activator or inhibitor), which is thresholded to yield a binary output.

Somewhat contrary to the previously mentioned studies, in [START_REF] Buchler | On schemes of combinatorial transcription logic[END_REF] mechanistic models of transcriptional control show that virtually any regulatory logic function can be implemented by cis-regulatory constructs, including e.g. XOR gates. Being neither canalyzing nor monotonous, these logic functions were deemed biologically meaningless in the studies mentioned above. In short, there is no generally accepted definition of "biologically meaningful".

Naturally, this matter becomes even more difficult, if nodes assume more than two discrete states. For this reason, we did not attempt to define the subset of biologically meaningful update rules. Rather, we delineated a characteristic feature -the single-switch condition -which is shared by (most) biologically meaningful update rules. We did so by analyzing manually curated multistate logical models as well as by studying a thermodynamic model of gene regulation. We do not claim that all single-switch functions are implemented in real gene regulatory networks. In other words, we deem our characterization of meaningful update rules rather conservative in the sense that the class of single-switch functions is likely too large, i.e. still contains biologically implausible rules. (Naturally, we cannot fully exclude the possibility that single update rules of real gene regulatory networks do not satisfy the single-switch condition, either.) This, of course, raises the question of how to define and impose further constraints.

A thermodynamic model of gene regulation

From a more general point of view, the presented thermodynamic model provides the biophysical foundation for the relation between discrete and continuous models of gene regulatory networks. Our modeling approach can be nicely related to previous studies:

• In the mathematical framework for gene regulatory networks described in [START_REF] Mestl | A mathematical framework for describing and analysing gene regulatory networks[END_REF] the probabilities π ω from (15) are called indicator functions; however, no derivation of these functions is presented.

• In the Boolean case, we can think of the discretized transcription rates r ω ∈ {0, 1}, ω ∈ {0, 1} K , as a K-variate Boolean function and the multilinear interpolation B I (rω) from ( 16) is equal to the Zhegalkin representation of this Boolean function [START_REF] Faisal | Polynomial models of gene dynamics[END_REF].

• The functions B P W (rω) from ( 19) give rise to an intermediate between discrete and continuous models: systems of piecewise linear differential equations as studied e.g. in [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF]; de [START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF].

• The established connection between continuous and discrete models of gene regulatory networks has already been used for the standardized continuous extension of Boolean models (Wittmann et al., 2009a,b).

The main reason for the complexity of transcriptional gene regulation are cooperative effects between transcription factors. Cooperativity exists, if a bound protein increases (positive cooperativity) or decreases (negative cooperativity) the affinity of the promotor for further proteins of the same kind (homotypic cooperativity) or of a different kind (heterotypic cooperativity).

In the derivation of formula (17) we made two simplifying assumptions regarding cooperative effects: We assumed strong homotypic cooperativity and neglected heterotypic cooperativity. Strong homotypic cooperativity allows us to assume a quasi-simultaneous binding of all ligands of the same kind. Intuitively speaking, we argue that once a binding site is occupied strong interactions between the ligands lead to a fast occupation of the remaining sites. Consequently, we can discard the configurations in which the binding sites are only partially occupied. This, in turn, greatly enhances the amenability of the resulting formulas -an important issue in 1910, when Hill conducted his famous studies [START_REF] Hill | The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves[END_REF].

Although nowadays computational intractability is no longer an issue, the Hill equation is still widely used due to its simplicity and variability. In fact, one can account for sequential ligand binding by allowing non-integer exponents. For example, in the case of hemoglobin, in which four oxygen molecules are known to bind with a high degree of positive cooperativity, the measured Hill coefficient ranges from 1.7 to 3.2 rather than 4 [START_REF] Hill | The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves[END_REF]. For a comprehensive review of the limitations of the Hill equation and possible more realistic extensions, see [START_REF] Weiss | The Hill equation revisited: uses and misuses[END_REF]. The multivariate case of these extensions, however, is still an open issue.

With respect to our second assumption, that of negligible heterotypic cooperativity, we remark that, for example, bioinformatics analyses of the Drosophila genome indeed revealed significant short-range homotypic clustering of binding sites but no systematic heterotypic clustering between binding sites of different factors [START_REF] Segal | Predicting expression patterns from regulatory sequence in Drosophila segmentation[END_REF]. Still, heterotypic cooperativity exists and we should be wary of this assumption. Luckily, we can argue that taking into account heterotypic cooperativity does not change the switch-like character of the overall transcription rate r. To see this, we fix all the regulator concentrations but one, w.l.o.g. [x 1 ]. Then we can write the overall transcription rate r in terms of

[x 1 ] r = c 1 + c 2 [x 1 ] n 1 c 3 + c 4 [x 1 ] n 1 ,
with constants c 1 , c 2 , c 3 , c 4 . This still is a sigmoidal interpolation between the transcription rate r = c 1 /c 3 for [x 1 ] = 0 and the transcription rate r ≈ c 2 /c 4 for large [x 1 ]. However, we can no longer express r in a form similar to (17), and also the determination of the switching thresholds becomes more involved.

Outlook

In continuation of our work, several further unrealistic features of KNs could be addressed. First, the synchronous updating according to an external clock is an over-simplification of biological reality. Note that we can think of (23) as a simple Euler approximation of (24).

The dynamical behavior of the continuous model has been shown to significantly differ from that of the Boolean network [START_REF] Glass | The logical analysis of continuous, nonlinear biochemical control networks[END_REF][START_REF] Mestl | A mathematical framework for describing and analysing gene regulatory networks[END_REF][START_REF] Edwards | Analysis of continuous-time switching networks[END_REF][START_REF] De Jong | Qualitative simulation of genetic regulatory networks using piecewise-linear models[END_REF][START_REF] Mochizuki | An analytical study of the number of steady states in gene regulatory networks[END_REF]. In particular, the detection and analysis of chaos in systems of piecewise linear differential equations is still largely an open issue, although it was addressed in several, mostly numerical, studies [START_REF] Mestl | Common chaos in arbitrarily complex feedback networks[END_REF][START_REF] Glass | Ordered and disordered dynamics in random networks[END_REF][START_REF] Kappler | Dynamics in high-dimensional model gene networks[END_REF].

Second, one could investigate the effect of different time-scales in the above system of piecewise linear differential equations ( 24) by multiplying the right-hand sides by time-scale parameters τ i . Several studies have shown that the integration of different time-scales increases the robustness of biological systems [START_REF] Gorban | Dynamical robustness of biological networks with hierarchical distribution of time scales[END_REF][START_REF] Rojdestvenski | Robustness and time-scale hierarchy in biological systems[END_REF]

= (p 1 + p 2 ) /2 it holds that 2p (1 -p ) = p 1 + p 2 - (p 1 + p 2 ) 2 2 > p 1 + p 2 - (p 1 + p 2 ) 2 2 - (p 1 -p 2 ) 2 2 = p 1 (1 -p 1 ) + p 2 (1 -p 2 ) .
Hence, the distribution

P S i (s) = ⎧ ⎨ ⎩ P S i (s) s / ∈ {s 1 , s 2 } p s ∈ {s 1 , s 2 }
yields a strictly larger value for p S i .

Appendix B.

Observe that for o > 0 the function O(o) is convex and monotonously increasing. To prove this, it suffices to show that o K + 1o K (1p S ) is convex and monotonously increasing. This can easily be verified by computing the first and second derivatives, In the latter case we have an additional fixed point, which attracts [0, 1). 4) for a frozen (blue) and a chaotic (green) KN.

Ko K-1 -(1 -p S ) Ko K-1 = Kp S o K-1 ≥ 0 and K(K -1)p S o K-2 ≥ 0,
The linear stability analysis of (4) about o yields:

dO do o = Kmax K=1 P in (K) Smax S=2 P nos (S) [K -K (1 -p S )] = Kmax K=1 P in (K) Smax S=2 P nos (S) Kp S = Kmax K=1 P in (K) Kp = pK ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 1 ordered regime = 1 critical boundary > 1 chaotic regime . Appendix C.
We substitute ( 14) in ( 5):

K 2w(1 -w) + (1 -w) 2 (1 -μ -1 ) ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 1 ordered regime = 1 critical boundary > 1 chaotic regime .
This is equivalent to

-K (1 + μ -1 ) w 2 + 2Kμ -1 w + K (1 -μ -1 ) -1 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 0 ordered regime = 0 critical boundary > 0 chaotic regime . (C.1)
The discriminant of this quadratic equation is given by

D = 4K K -(1 + μ -1 ) and D ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 0 if K < 1 + μ -1 = 0 if K = 1 + μ -1 > 0 if K > 1 + μ -1 .
From (C.1) it follows that for K < 1+μ -1 the system is always frozen. To see this, note that -K (1 + μ -1 ) < 0. For K ≥ 1 + μ -1 the quadratic equation from (C.1) has two (possibly coinciding) solutions

w crit 1 2 K, μ -1 = 1 μ -1 + 1 μ -1 ∓ 1 - μ -1 + 1 K First, observe that for μ -1 > 0, w crit 2
< 1 is always a valid solution. Moreover, we have

w crit 1 > 0 ⇔ 1 - μ -1 + 1 K < μ -1 ⇔ K < 1 1 -μ -1 ,
where in the case μ -1 = 1, "K < ∞" is true for any K. Observe that for

μ -1 < 1, 1 1 -μ -1 > 1 + μ -1 ,
so there always is a range for K, in which (C.1) has two distinct solutions in (0, 1). Once more considering that in (C.1) -K (1 + μ -1 ) < 0, we finally obtain

-K (1 + μ -1 ) w 2 +2Kμ -1 w+K (1 -μ -1 )-1 ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ < 0 , w < w crit 1 or w > w crit 2 = 0 , w = w crit 1 or w = w crit 2 > 0 , w crit 1 < w < w crit 2 . Appendix D.
Using μ -1 = 1/(S -1) we can rewrite the classification from section 3.1 in terms of S. We also write w crit (K, S) instead of w crit (K, μ -1 ).

(a) For K < S/(S -1) the system is always frozen.

(b) For K = S/(S -1) the system is critical if and only if

w = w crit S S -1 , S = 1 S ,
i.e. P S is the uniform distribution, otherwise it is frozen.

(c) For S = 2 and K > S/(S -1) the critical boundary is described by

w crit 1 2 K, 2 = 1 2 1 ∓ 1 - 2 K ∈ (0, 1) . (D.1) For w crit 1 K, 2 < w < w crit 2 K, 2
the system is chaotic, otherwise it is frozen.

(microstates); we describe them by Boolean vectors ω of length K k=1 n k , the b-th component ω(b) indicating if the respective binding site is occupied. (The reader is kindly asked to ignore that this definition of ω conflicts with the one given in the main text, we will resolve this discrepancy below.)

For each microstate ω we denote the number of binding sites occupied by regulator x k by #ω(k).

Each configuration ω has a certain Gibbs free energy G ω and the probability of a certain microstate ω 0 is given by

π ω 0 = exp (-ΔG ω 0 /RT ) K k=1 [x k ] #ω 0 (k) ω exp (-ΔG ω /RT ) K k=1 [x k ] #ω(k) , (E.1)
where Δ signifies energy differences with respect to the state (0, 0, . . . , 0), R is the gas constant, T the temperature, and brackets indicate concentrations of free ligand. Since we may assume the number of ligands to be larger than the number of promotors by several orders of magnitude, we can neglect the difference between concentration of free ligand and total ligand concentration.

Formula (E.1) can also be written in terms of the association constants

A ω , which are related to the free energies via ΔG ω = -RT ln A ω . Observe that ΔG (0,0,...,0) = 0 implies A (0,0,...,0) = 1. We obtain k) . For simplification, one commonly assumes that either no protein or n 1 proteins are bound to the promotor, i.e. one discards all microstates except 0 = (0, 0, . . . , 0) and 1 = (1, 1, . . . , 1). A possible justification of this assumption is a strong positive homotypic cooperativity. Indeed, according to (E.1)

π ω 0 = A ω 0 K k=1 [x k ] #ω 0 (k) ω A ω K k=1 [x k ] #ω(
and (E.3) we can neglect all π ω except π 0 and π 1 if -ΔG hom 1 is sufficiently large. This issue is further discussed in section 6.4.

Equation (E.2) then becomes

π 0 = 1 1 + A 1 [x 1 ] n 1 , π 1 = A 1 [x 1 ] n 1 1 + A 1 [x 1 ] n 1 ,
or, in terms of the dissociation constant D 1 = 1/A 1 ,

π 0 = D 1 D 1 + [x 1 ] n 1 , π 1 = [x 1 ] n 1 D 1 + [x 1 ] n 1 .
These are the well-known negative and positive Hill functions [START_REF] Hill | The possible effects of the aggregation of the molecules of haemoglobin on its dissociation curves[END_REF],

π ω 0 = A ω 0 k|ω 0 (k)=1 [x k ] n k ω A ω k|ω(k)=1 [x k ] n k .
(E.4)

In general, it holds that

ΔG ω = k|ω(k)=1 ΔG e k + ΔG het ω ,
where ΔG e k is the free energy of binding the k-th regulator polymer and ΔG het ω is the free energy of interactions between the different regulators (heterotypic cooperativity). In the following, we neglect heterotypic cooperativity, all ΔG het ω = 0, for a discussion see section 6.4. In terms of the association constants this implies

A ω = k|ω(k)=1
A e k .

(E.5)

Hence, the denominator in (E.4) factorizes

ω A ω k|ω(k)=1 [x k ] n k = ω k|ω(k)=1 A e k [x k ] n k = K k=1 (1 + A e k [x k ] n k ) . (E.6)
Using (E.5) and (E.6) we can rewrite (E.4) as

π ω 0 = k|ω 0 (k)=1 A e k [x k ] n k K k=1 (1 + A e k [x k ] n k ) = k|ω 0 (k)=1 A e k [x k ] n k 1 + A e k [x k ] n k • k|ω 0 (k)=0 1 1 + A e k [x k ] n k = k|ω 0 (k)=1 h + ([x k ] , n k , A e k ) • k|ω 0 (k)=0 h -([x k ] , n k , A e k ) .
Replacing A e k by A k and dropping the subscript 0 we obtain (15).

Figure 1 :

 1 Figure 1: (A) The case of unbiased update rules. Critical connectivities K crit `S´f rom (11) (blue 'x') and K crit (Smax) from (12) (green '+'). (B) The case of biased update rules. Critical boundaries w crit (K, S) from section 3.2. For S = 2 the two branches w crit 1/2 `K, 2 ´from (D.1) are shown. For S = 3 the critical boundary consists of three parts: w crit 1/2 `K, 3 ´from (D.2) and w crit `K, 3 ´from (D.3). (C) Plots of the overlap o(t) from (3) between two initially differing time courses for a frozen, critical and chaotic network. Each network has N = 100 nodes, each node has S = 3 states and is connected to K = 4 inputs. The update rules are constructed according to distribution (13). For the frozen/critical/chaotic network we set w = 0.95 / w = w crit (4, 3) ≈ 0.8604 from (D.3) / w = 1/3. The initial overlap was chosen to be o(0) = 0.9. (D,E,F) The first 50 time steps of evaluations of the frozen, critical and chaotic networks from (C); only the first 25 nodes are shown. The three states are color coded in black, grey and white.Note that subfigures C,D,E,F all have the same X-axis, namely time running from 1 to 50 as indicated in (F).

  Figure 1F. (We know, of course, that ultimately also the chaotic network will fall into a limit cycle due to the finite size of the state space.) If we discard the first 10 time steps, the overlap in the critical network oscillates around 0.88 and fluctuates around 0.34 in the chaotic network. The latter ratio agrees well with the expected overlap of two randomly drawn configurations.

Figure

  Figure 2A shows an example of a random multistate update rule. Is such a function really plausible from a biological point of view? To answer

Figure 2 :

 2 Figure 2: A Random update rule depending on two inputs, x 1 ∈ {1, 2, 3, 4} and x 2 ∈ {1, 2, 3}. The function values were drawn randomly from {1, 2, 3, 4, 5}. B Single-switch function depending on two inputs, x 1 ∈ {1, 2, 3, 4} and x 2 ∈ {1, 2, 3}. The function possesses one switching point along the x 1dimension (between 2 and 3) and one switching point along the x 2 -dimension (between 1 and 2). C Positive Hill function and its approximation by a Heaviside step function. D,E,F Derivation of the single-switch function shown in (B) from transcription rates r (0,0) = 0, r (1,0) = 0.25, r (0,1) = 0.75, r (1,1) = 1 as discussed in section 4.4. From left to right: B I (rω ) from (16), overall transcription rate r from (17) and B P W (rω ) from (19). G Schematic representation of transcriptional gene regulation. Clusters of identical binding sites allow for the homotypicly cooperative binding of the transcription factors x 1 , x 2 , . . . , x K to the promotor region of gene x. Heterotypic cooperativity is not taken into account. H Critical connectivities K crit for random and single-switch update functions in terms of a fixed number of states S. Blue markers ('x') show the critical connectivity from (11) for fully random update rules; green markers ('+') show the critical connectivity from (22) for single-switch functions.

  are positive and negative Hill functions, respectively. Brackets indicate ligand concentrations, and A k is the association constant describing the bonding affinity between the k-th regulator and its binding site.We denote the transcription rate in each promotor configuration ω by r ω and let B I (rω) be the K-dimensional multilinear interpolation of the points (

..

  n 1 , n 2 ) we obtain the overall transcription rate r from (17) shown in Figure 2E. Finally, replacing the Hill functions by Heaviside step functions yields the piecewise linear function B P W (rω) from (19) shown in Figure 2F. Now, let us say that x 1 regulates two other genes, one at a higher and one 22 at a lower threshold than D micro 1 In this case, we would allow 4 discrete states for x 1 separated by these three thresholds and D micro 1 would lie between the second and the third discrete state. Similarly, let x 2 regulate one other gene at a higher threshold than D micro 2 This would require three discrete states for x 2 and D micro 2

  respectively. Moreover, note that p S < 1 implies O(0) > 0. Hence, there are two possibilities: Either the fixed point o is (globally) attractive and no further fixed points exist, or o is repellent and O has exactly one additional (attractive) fixed point. These two possibilities are visualized in Figure B.3. We observe that for the blue curve the fixed point o is attractive (even globally), while for the green curve it is repellent.

  Figure B.3: Schematic plots of O(o) from (4) for a frozen (blue) and a chaotic (green) KN.

  by considering the case K = 1. We denote the free energy of binding a protein to the b-th binding site by ΔG e b , where e b is the b-th unit vector of length n 1 . The free energy ΔG ω of a microstate ω can then be written as ΔG ω = b|ω(b)=1 ΔG e b + ΔG hom ω , (E.3) i.e. as the sum of the free energies ΔG e b of the single protein-DNA bindings and the free energy ΔG hom ω of interactions between the proteins (homotypic cooperativity). Clearly, ΔG hom e b = 0 for all unit vectors e b .

  which are often used to model transcriptional gene regulation. Typically, they are stated in terms of the microscopic dissociation constant D micro 1 = n 1 D 1 , which indicates the concentration level for half-maximal activation. The basic idea behind discrete models of gene regulatory networks is to approximate the sigmoidal Hill functions by Heaviside step functions, cf. Figure 2C. Let us now consider the case of general K. For simplification, we again assume that each transcription factor x k , 1 ≤ k ≤ K, can bind only as a polymer consisting of n k monomers. This limits the number of accessible microstates to 2 K ; by abuse of notation we denote them again by vectors ω ∈ {0, 1} K , the k-th component ω(k) indicating whether or not the k-th regulator polymer is bound. With this new definition of ω, which now agrees with the one given in the main text, equation (E.2) reads

  . So far, however, authors almost exclusively treated the special case of equal timescale parameters. The investigation of the general system will require the development of new analytical and numerical approaches. Still, we believe that the effort is worthwhile. Studying the effect of different time-scales on the self-organizing property of random networks will be interesting from both, a biological as well as a general systems-theoretic point of view.Third, it remains to be investigated whether results about piecewise linear differential equations remain valid, if we do not approximate the sigmoidal functions in (17) by step functions. At least for sufficiently steep sigmoidals we expect this to be the case. However, from a biological point of view, the steepness is limited by the number of binding sites. This raises the question

	of a global maximum follows from the extreme value theorem. Now assume
	that p 1	def = P S i (s 1 ) = p 2	def = P S i (s 2 ). Then, for p	def
	to what extent the "mathematical lower bound" and the "biological upper
	bound" can be reconciled.

Acknowledgements

The authors wish to thank Florian Blöchl (Helmholtz Center, Munich) for critical reading of the manuscript and useful discussions. Funding by the Helmholtz Alliance on Systems Biology (project CoReNe) is gratefully acknowledged.

Appendix A.

We proof that the expression for p S i from equation (1) has a unique global maximum if P S i is the discrete uniform distribution from (2). The existence (d) For S > 2 and (S -1)/(S -2) > K > S/(S -1) the critical boundary is described by

the system is chaotic, otherwise it is frozen. For K = (S -1)/(S -2), w crit 1 (K, S) = 0 is not a valid solution for w as we require w ∈ (0, 1). Note that for w = 0, P S from (13) would essentially become a uniform distribution on 1, 2, . . . , S -1. We could replace S by S = S -1 and choose a uniform distribution for P S . This new system falls into category (b) and we find that it is in fact critical.

(e) For S > 2 and K ≥ (S -1)/(S -2) the critical boundary is described by [START_REF] Solé | Phase transitions in random networks with multiple states[END_REF], for larger values of w the system is frozen, otherwise it is chaotic.

Appendix E.

Let us assume a situation as outlined at the beginning of section 4.2. We begin our computations without any assumptions regarding cooperative effects. Hence, there are 2

Let us consider the i-th node, 1 ≤ i ≤ N , and its k-th input, 1 ≤ k ≤ K i .

The probability q S ik that for two different values x ik = x ik of this input we have h step (x ik , θ) = h step (x ik , θ), where θ ∈ Σ ik \ {0} is a randomly chosen threshold, is given by

To see this, think of x ik and x ik as i.i.d. random variables distributed according to P S ik . The bracketed expression is the product of the probability for "x ik = s" and the conditional probability for "x ik ≥ θ" given "x ik = s".

The sum over s thus gives the probability for "x ik < θ and x ik ≥ θ" given "x ik = x ik ". Twice this sum is the probability that x ik and x ik lie on different sides of θ and thus give different values when plugged into h step (•, θ). Computing the weighted sum over all non-degenerate thresholds θ finally yields the above formula (F.1).

In the situation of section 5.2 formula (F.1) simplifies to

It is interesting to compare this to the continuous situation, where the two input values and the threshold are drawn from some continuous interval, w.l.o.g. the unit interval [0, 1]. Intuitively one might suspect that the probability q cont for the two input values to lie on different sides of the threshold is 1/2. However, q cont computes to

and agrees with the limit of q from (F.2) as S → ∞.


Gabriel Sargent
email: gabriel.sargent@irisa.fr

Frédéric Bimbot
email: frederic.bimbot@irisa.fr

Emmanuel Vincent
email: emmanuel.vincent@inria.fr

A REGULARITY-CONSTRAINED VITERBI ALGORITHM AND ITS APPLICATION TO THE STRUCTURAL SEGMENTATION OF SONGS

This paper presents a general approach for the structural segmentation of songs. It is formalized as a cost optimization problem that combines properties of the musical content and prior regularity assumption on the segment length. A versatile implementation of this approach is proposed by means of a Viterbi algorithm, and the design of the costs are discussed. We then present two systems derived from this approach, based on acoustic and symbolic features respectively. The advantages of the regularity constraint are evaluated on a database of 100 popular songs by showing a significant improvement of the segmentation performance in terms of F-measure.

INTRODUCTION

Music structure is one of the properties which contributes to the characterization of a music piece. It describes its temporal organization at a high level, by means of segments labeled according to their musical content and their relationships with one another. The automatic structural segmentation of songs is generally addressed by analyzing the homogeneity and the repetitiveness of the musical content over time (timbre, harmony, rhythm, melody).

Recent work [START_REF] Bimbot | Decomposition into autonomous and comparable blocks : a structural description of music pieces[END_REF] proposes a single-level definition of the structure of a music piece based on a regularity assumption. It implies the prevalence of one (or a few) typical segment duration(s) within each song, i.e. structural pulsation period(s). Indeed, a large part of western popular music is built on musical patterns (rhythmic cells, chord progressions, melodies...) which show cyclic behaviors and which are fully or partly repeated over time. This induces some sort of regularity in the structure of songs.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. c 2011 International Society for Music Information Retrieval.

The present work is based on this regularity assumption in music. We introduce a general segmentation framework, which consists of an optimization method to find the best segmentation combining the similarities and/or the contrasts in musical content and the regularity of the segments. An implementation of this method is proposed by means of a Viterbi algorithm.

A similar segmental Viterbi algorithm was briefly sketched in [START_REF] Levy | New methods in structural segmentation of musical audio[END_REF] in the context of a probabilistic model (segmental HMM). In this paper, we make it more explicit and we extend it to any type of cost function. This makes it possible to exploit combinations of clustering-based and similarity-matrix-based approaches and to a wider variety of situations outside the probabilistic framework. We also discuss the importance of the regularity cost in the estimation of the segment boundaries, and provide experimental results with several choices for the two terms of the segmentation cost.

The structure of the paper is as follows. In section 2 we present the general music segmentation method, without considering a particular musical feature or temporal scale. Section 3 describes its implementation by means of a Viterbi algorithm, and discusses the expression of segmentation costs. In section 4, after briefly reviewing former work on music structure, we apply the proposed segmentation method to this particular problem. We then present two structural segmentation systems based on the algorithm developed above. Section 5 evaluates the effect of the incorporation of regularity constraints thanks to the evaluation of these systems on the RWC popular music database [START_REF] Goto | RWC Music Data-base : RWC music database : Popular, Classical, and Jazz Music Databases[END_REF].

GENERAL APPROACH

This section presents a general method for the temporal segmentation of music pieces, when regularity assumptions can be hypothesized on the segment length. It consists of an optimization process where the optimal segmentation is searched simultaneously considering the properties of the data and the regularity of the segmentation.

A music piece X can be described as a sequence of N features {x t } 1≤t≤N along a particular temporal scale (e.g. frames, or beats...). We denote X tj ti = {x t } ti≤t<tj the sequence of features associated to the temporal interval [t i , t j [.

Let us define a segmentation S = {s k } 1≤k≤n of X as a sequence of n intervals s k = [t k , t k+1 [, with the following conventions :

• t 1 = 1 < ... < t k < ... < t n < t n+1 = N + 1, • s 0 = [t 0 , t 1 [= [0, 1[, for the algorithm initialization, • m k = t k+1 -t k is the length of s k .
We aim at finding the optimal segmentation, by minimizing a certain cost function.

We assume that the cost function C can be written as

C(S) = n k=1 Γ(s k ) (1)
with

Γ(s k ) = Φ(s k ) + λ(τ )Ψ(s k ) + ǫ (2)
where

• Φ(s k
) is a content-based segmentation cost, which takes low values when the sequence of features in s k is likely to correspond to a structural segment. This cost can be described according to different families of functions, like change detection functions or similarity functions. It can also, for instance, be derived from a probabilistic function P (s k ), as -logP (s k ).

• Ψ(s k ) is a regularity cost. We consider that the regularity of a segmentation depends on the deviation of the length of its segments to a prior reference length τ called the structural pulsation period (as a consequence, Ψ(s k ) decreases as m k approaches τ ). Note that, if the values of m k are expected to follow a particular distribution π(m k ) around τ , Ψ(s k ) can be set as Ψ(s k ) =log(π(m k )).

• λ(τ ) is a balance parameter between these two costs.

• In practice, we add a small constant ǫ > 0 to give a slight advantage to longer segments in the case where Φ and Ψ would be equivalent for several segmentations.

IMPLEMENTATION

This section presents an implementation of the approach presented above, and describes possible choices of cost functions Φ, Ψ and parameter λ. Let s t,h be the interval corresponding to X t t-h = [x t-h , x t [, the set of features which precede the temporal index t within a window of length h. We denote H as the maximal window length considered1 .

• Initialization (t = 1)

We set S 1 = {[0, 1[} and C 1 = 0. • For t = 1 : N -1 We consider {s t,h } 1≤h≤H ′ , with H ′ = min(t -1, H)
as the set of admissible predecessors for temporal index t.

The optimal segmentations {S t-h } 1≤h≤H ′ ending at indexes {t -h} 1≤h≤H ′ are assumed to be known, as well as their associated cumulative costs

{C t-h } 1≤h≤H ′ .
Then, the best partial segmentation S t is built by choosing the extension of the former partial segmentation S t-h with the lowest cost. We evaluate respectively :

1. Γ(s t,h ) for 1 ≤ h ≤ H ′ , 2. b(t) = argmin 1≤h≤H ′ {C t-h + Γ(s t,h )}, 3. C t = C t-b(t) + Γ(s t,b(t) ) We can note that S t = S t-b(t) ∪ {S t,b(t) }.
The optimal segmentation for X, noted S opt with cost C N +1 , is obtained by backtracking the optimal predecessors stored in b(t). The associated temporal indexes {t k } 1≤k≤nopt are then found thanks to the following recursion :

1. t nopt+1 = b(N + 1), 2. t k = b(t k+1 ), for 1 ≤ k ≤ n opt .
n opt is the number of boundaries of S opt , obtained after this backtracking process.

Design of the cost functions

Content-based segmentation cost Φ

The objective of the content-based segmentation cost is to evaluate a set of segments according to the redundancy of their content. Segmentations with lower costs are expected to consist of segments built on the same musical patterns. Different families of functions can be considered, like abrupt change detection criteria or similarity functions.

Abrupt change detection criteria assign a low cost to segments associated to probable boundaries. In automatic structure inference, [START_REF] Cooper | Media segmentation using selfsimilarity decomposition[END_REF] uses for example a "novelty function" based on the analysis of the local homogeneity of the song over time.

Similarity functions aim to assign a low cost to segments made of sequences of features repeated elsewhere in the song. We can define such a function as

Φ(s k ) = min θ∈Z k {φ(X t k +m k t k , X θ+m k θ )}. (3)
The lowest dissimilarity φ is taken between the sequence of features X t k +m k t k from s k (of length m k ) and any other sequence of the same length contained in a portion Z k of X. In particular, Φ(s k ) = 0 when the sequence of features of s k is exactly repeated elsewhere in

Z k , Φ(s k ) > 0 otherwise. Z k = [1, t k-m k ] ∪ [t k+m k , N
] can be chosen to avoid intra-segment comparisons. In the case of a binary dissimilarity, where a song is described as a sequence of symbolic features, the following function can be chosen :

φ(X t k +m k t k , X θ+m k θ ) = m k -1 p=0 1 -δ(x t k +p , x θ+p ), ( 4
)
where δ is Kronecker's delta (equals 1 when arguments have the same value, 0 otherwise). More generally any nonbinary function can be used in equation (3).

Regularity cost Ψ

The regularity cost Ψ of a segmentation is based on the measure of the deviation between the length of its segments from a reference length τ . It can show the following properties :

1. Ψ(τ ) = 0, 2. Ψ(m k ) > 0, taking higher values as the segment length m k moves away from τ .

A lot of functions can satisfy these properties. We consider two categories of functions : convex and non-convex functions. As non-convex functions verify the property :

Ψ(τ ) + Ψ(τ + ∆) < Ψ(τ + ∆ 1 ) + Ψ(τ + ∆ 2 ) (5)
with As an illustration, we consider the following family of symmetric functions derived from the l α norm :

∆ = ∆ 1 + ∆ 2 , (6)
Ψ α (m k ) = | m k τ -1| α ( 8
)
m k is the length of interval s k , and α controls the convexity of the function (we have a non-convex function if 0 < α < 1, and a convex one if α ≥ 1). Figure 2 shows Ψ α for α = {0.5, 1, 2}.

Balance parameter λ

We consider that λ depends on τ as the probability of having irregular segments grows with the number of segments, and therefore with the inverse of τ . We choose the linear relation λ(τ ) = λτ , where λ is a constant parameter to be tuned.

APPLICATION TO THE STRUCTURAL SEGMENTATION OF SONGS

The work presented in section 3 is primarily intended to the structural segmentation of songs. Automatic music structure inference is a difficult task, because the problem to be solved is usually ill-posed. Moreover, it requires the analysis and the complex combination of features and criteria through the development of sophisticated metrics and algorithms. In this section, we review briefly the main state-of-the-art methods for automatic structural segmentation of songs, before describing two structural segmentation systems implemented from the proposed method.

State-of-the-art

Different approaches have been proposed to the problem of automatic structure inference. They generally use acoustic features, such as Mel-Frequency Cepstral Coefficients (MFCCs) and Chroma vectors, which characterize the instrumental timbre and the harmonic content respectively. Other features are described in [START_REF] Turnbull | A supervised approach for detecting boundaries in music using difference features and boosting[END_REF], [START_REF] Barrington | Modeling music as a dynamic texture[END_REF], and [START_REF] Jensen | Multiple scale music segmentation using rhythm, timbre and harmony[END_REF]. Structural segments are assumed to show stable instrumentation (often associated to homogeneous timbre) and therefore to appear as blocks with specific textures in similarity matrices [START_REF] Cooper | Media segmentation using selfsimilarity decomposition[END_REF][START_REF] Paulus | Music structure analysis by finding repeated parts[END_REF], or sequences of similar states in Hidden Markov Models (HMMs) [START_REF] Logan | Music summarization using key phrases[END_REF]. Repeated harmonic progressions can be detected by localizing the sequences of high similarity coefficients in subdiagonals of the chroma-based similarity matrix [START_REF] Goto | A chorus-section detecting method for musical audio signals[END_REF]. Other approaches, like HMMs [START_REF] Levy | New methods in structural segmentation of musical audio[END_REF][START_REF] Peeters | Toward automatic music audio summary generation from signal analysis[END_REF], or more recently Non negative Matrix Factorization [START_REF] Weiss | Identifying Repeated Patterns in Music Using Sparse Convolutive Non-Negative Matrix Factorization[END_REF] are also used for the recognition of repeated harmonic patterns. Some methods use dynamic programming : Shiu et al. interpret the chroma-based similarity matrix as a time-state representation and use the Viterbi algorithm to find the path with the highest score in terms of similarity through it [START_REF] Shiu | Similarity matrix processing for music structure analysis[END_REF]. A constraint is set to give priority to the diagonal direction for the path, and implicitly influence the length of the estimated structural segments.

Some other approaches combine these content-based methods by means of optimization problems, as in [START_REF] Jensen | Multiple scale music segmentation using rhythm, timbre and harmony[END_REF][START_REF] Paulus | Music structure analysis by finding repeated parts[END_REF]. A more detailed state of the art is available in [START_REF] Paulus | Audio-based music structure analysis[END_REF].

In the following section, we present two systems that infer the structural segmentation of a song, incorporating the idea of "structural pulsation period"2 .

Presentation of the systems

These systems perform a structural segmentation of songs combining content-based segmentation under a regularity constraint by means of the Viterbi algorithm presented in section 3.1. System 1 uses acoustic features to compute change detection criteria and estimates the main structural pulsation period τ from the audio. System 2 analyzes symbolic features, uses a similarity function and prior knowledge of τ (fixed at 32 beats). As features are considered at the beat scale, a beat detection system is needed. We evaluate for these 2 systems the impact of incorporating a regularity constraint on the relative performance of the segmentation.

System 1 : combination of change detection criteria on acoustic features

The system we consider is the one described in [START_REF] Sargent | A structural segmentation of songs using generalized likelihood ratio under regularity assumptions[END_REF]. In this paper, we consider variants of this system both with and without the regularity constraint in order to analyze its impact on structural segmentation inference. The contentbased segmentation cost is based on 3 statistical criteria which measure for each temporal index the likelihood ratio of a structural segment boundary. This criterion combines instrumental changes, short events and contrastive patterns over time.

The criteria are combined in a weighted sum to form what we name here the content-based segmentation cost. A linear regularity cost function is used to perform the Viterbi approach described in section 3.1, to find the segmentation with lowest cost. The main structural pulsation period of the song is estimated by a Fourier transform on the instrumental change criterion.

System 2 : similarity function on symbolic features

It is interesting to consider symbolic features for structure inference as other means of music description. The joint use of various features in a global and versatile retrieval system may increase the accuracy of the estimated segmentation [START_REF] Vincent | A roadmap towards versatile MIR[END_REF]. The symbols can be obtained for instance from a score of the piece. System 2 uses chords estimations to compute the similarity function described with the equations ( 3) and ( 4) of section 3.2.1. Each chord class is associated to a different symbol, to obtain a quite neutral symbolic description of the song. The size of the alphabet of symbols we use is the number of chord classes used by the chord estimator (e.g. 24 classes for major and minor chords). Each symbol corresponds to a duration of 2 beats, in order to be consistent with the temporal scale used in [START_REF] Bimbot | Decomposition into autonomous and comparable blocks : a structural description of music pieces[END_REF].

The structural pulsation period value τ is considered as prior knowledge and used in the regularity cost Ψ α of equation (8), section 3.2.2. The content-based cost and the regularity cost are then combined using equations ( 1) and ( 2) from section 2, and the segmentation with lower cost is found using Viterbi algorithm from section 3.1.

EVALUATION

Evaluation database

The algorithms have been evaluated using the RWC popular music database [START_REF] Goto | RWC Music Data-base : RWC music database : Popular, Classical, and Jazz Music Databases[END_REF], and the set of reference annotations provided by [START_REF] Bimbot | Decomposition into autonomous and comparable blocks : a structural description of music pieces[END_REF], used in MIREX 2010. This database consists of 100 songs written and produced for research purposes.

Evaluation metrics

The evaluation of the segmentation is done by Precision (P ), Recall (R) and F-measure (F ) metrics. Let s R be the set of reference boundaries (annotations) and s E the set of estimated ones, they are respectively defined as :

P = |s E ∩ s R | |s E | ; R = |s E ∩ s R | |s R | ; F = 2P R (P + R) . (9)
The matching of reference and estimated boundaries is performed within particular tolerance windows. We consider 0.5 s and 3 s as in MIREX 2010. Note that each boundary is used only once during the matching process.

Feature extraction and algorithm parametrization

System 1 (which uses change detection criteria) uses 20 MFCCs (including the 0th coefficient), extracted from Tolerance window = 0.5 sec. system α λ P (%) R(%) F (%) [START_REF] Ueda | HMM-based Approach for Automatic Chord Detection Using Refined Acoustic Features[END_REF], and uses the downbeat annotations available with the RWC database 5 . The reference annotations show that more than 80% of the songs have a main structural pulsation of 32 beats. We will then use τ = 32 beats as prior knowledge for our evaluation, and H = 3τ as the maximal number of admissible predecessors for each temporal unit.

A preliminary study on a subset of RWC popular was carried out to identify reasonable values of λ which fall within the interval [0, 1]. Three values of α are chosen to consider regularity costs functions with different convexities : a non-convex regularity cost function (α = 0.5), a convex regularity cost function (α = 2), and the intermediate case α = 1.

Results

The values gathered in Tables 1 and2 for System 1 and 2 show that the overall mean F-measures increase significantly when the regularity cost is introduced.

Figure 3 shows the average F-measure obtained with System 2 for the 3 regularity costs mentioned in 5.3. The values of λ corresponding to optimal performance appear in Table 2 for each case. The value of α has an impact on the accuracy of the estimated boundaries : it can be seen that, for a small tolerance, a non-convex regularity cost function gives better boundary accuracy than a convex one. This can be explained by the fact previously mentioned, that the convex case (α = 2) tends to spread structural irregularities (devi- ations from the ideal segmentation with segments of length τ ) across several structural segments. On the contrary, the non-convex case (α = 0.5) tends to concentrate them on a few segments. These results therefore show not only the advantage of the regularity constraint but also the importance of the fine properties of the corresponding cost function.

As a point of comparison, the best system in structural segmentation at MIREX 20106 (MND1) obtained Fmeasures of 35.9% and 60.5% (for tolerance windows of 0.5 s and 3 s respectively) on the same database. Note however that System 2 relies on a manual annotation of the downbeats.

CONCLUSION

The work presented in this paper has highlighted the relevance of incorporating a regularity constraint in the task of structural segmentation. Even with very basic cost functions as the ones considered in the present work, the very existence of the regularity constraint favors the retrieval of a well-defined solution. The Viterbi implementation, which we have detailed, allows a fast calculation of the optimal solution, and it can be applied in a generic way to any type of cost function.

The corresponding Matlab code will be made available to the MIR community 7 for enabling further experimental investigation within diverse structural segmentation systems and possibly for other tasks in MIR where the regularity constraint can be meaningful.

Figure 1 .

 1 Figure 1. Admissible predecessors for t and their costs

Figure 2 .

 2 Figure 2.Examples of regularity costs Ψ α for α = {0.5, 1, 2} and τ = 32

Figure 3 .

 3 Figure 3. Evolution of the average F-measures of System 2 on the database described in 5.1, as a function of balance parameter λ, for 3 types of regularity cost function (Ψ α={0.5,1,2} , τ = 32).

Table 1 .

 1 Average Precision (P ), Recall (R), and F-measure (F ) for two versions of System 1 on the RWC pop database. frames of length 23.2 ms, and a hop size of 11.6 ms (using scripts from MA toolbox by Beth Logan and Malcolm Slaney 3 ). Chroma vectors (12 coefficients) are extracted from frames of length 92.9 ms, and a hop size of 23.2 ms. Chroma vectors and beats estimation are computed thanks to LabRosa scripts 4 . System 2 (based on a similarity function) inputs the chords transcriptions obtained by the algorithm from Ueda et al., described in

	1	-	0	10.1	53.6	17.0
		1 0.5	23.9	24.3	23.8
		Tolerance window = 3 sec.
	system α	λ	P (%) R(%) F (%)
	1	-	0	16.8	89.3	28.2
		1 0.5	61.2	63.0	61.4

Table 2 .

 2 Average Precision (P ), Recall (R), and F-measure (F ) for System 2 (optimally tuned, considering λ ∈ [0 : 1]), on the database described in 5.1.

		Tolerance window = 0.5 sec.
	system	α	λ	P (%) R(%) F (%)
	2	-	0	17.9	31.9	22.0
		0.5 0.30	37.7	34.8	35.6
		1	0.30	34.7	32.3	33.0
		2	0.95	29.3	26.8	27.5
		Tolerance window = 3 sec.
	system	α	λ	P (%) R(%) F (%)
	2	-	0	36.1	64.7	44.5
		0.5 0.15	63.1	63.1	62.0
		1	0.15	63.4	64.1	62.7
		2	0.60	64.5	60.0	61.2

Typically, H = N , but smaller values can be used (e.g. multiples of τ ).

This can be seen as a way to constrain the ill-posed problem of structural segmentation towards a well-defined solution.

http://www.ofai.at/elias.pampalk/ma/documentation.html

http://labrosa.ee.columbia.edu/projects/coversongs/

http://staff.aist.go.jp/m.goto/RWC-MDB/

http://nema.lis.illinois.edu/nema out/mirex2010/results/struct/mirex10 /summary.html

http://www.irisa.fr/metiss/logiciel/

ACKNOWLEDGEMENTS

The authors would like to thank Yushi Ueda and Nobutaka Ono for their help in the collection of chord transcriptions used in this article. This work was partly supported by the Quaero project 8 funded by Oseo and by the associate team VERSAMUS 9 funded by INRIA.