Théotime Grohens
email: theotime.grohens@insa-lyon.fr

Sam Meyer

Guillaume Beslon

A Genome-Wide Evolutionary Simulation of the Transcription-Supercoiling Coupling

Keywords: DNA supercoiling, gene transcription, genome, evolution

DNA supercoiling, the level of under-or overwinding of the DNA polymer around itself, is widely recognized as an ancestral regulation mechanism of gene expression in bacteria. Higher levels of negative supercoiling facilitate the opening of the DNA double helix at gene promoters, and thereby increase gene transcription rates. Different levels of supercoiling have been measured in bacteria exposed to different environments, leading to the hypothesis that variations in supercoiling could be a response to changes in the environment. Moreover, DNA transcription has been shown to generate local variations in the supercoiling level, and therefore to impact the transcription rate of neighboring genes.

In this work, we study the coupled dynamics of DNA supercoiling and transcription at the genome scale. We implement a genome-wide model of gene expression based on the transcriptionsupercoiling coupling. We show that, in this model, a simple change in global DNA supercoiling is sufficient to trigger differentiated responses in gene expression levels via the transcription-1 supercoiling coupling. Then, studying our model in the light of evolution, we demonstrate that this non-linear response to different environments, mediated by the transcription-supercoiling coupling, can serve as the basis for the evolution of specialized phenotypes.

Introduction

The DNA molecule is a double-stranded polymer of nucleotides that plays a fundamental role in life.

It is shaped as a double helix which rotates around itself at a rate of around one turn per 10.5 base pairs [START_REF] Krogh | Impact of Chromosomal Architecture on the Function and Evolution of Bacterial Genomes[END_REF]. However, when subject to physical forces, it can become overwound or underwound, or writhe around itself, in a process known as DNA supercoiling; the supercoiling level σ is measured as the density of extra turns (or coils) per base pair [START_REF] Duprey | The regulation of DNA supercoiling across evolution[END_REF].

In bacterial cells, DNA is usually slightly underwound [START_REF] Lal | Genome scale patterns of supercoiling in a bacterial chromosome[END_REF], with a negative value of σ. The basal supercoiling level depends on the species, but a typical value is σ basal ≈ -0.066 in Escherichia coli [START_REF] Crozat | Long-Term Experimental Evolution in Escherichia coli . XII. DNA Topology as a Key Target of Selection[END_REF]. DNA supercoiling is tightly regulated by a class of enzymes called topoisomerases. The main topoisomerases are topoisomerase I and gyrase: gyrase uses ATP to maintain DNA in a negative supercoiling state by adding negative coils, while topoisomerase I relaxes supercoiling and does not need ATP [START_REF] Martis | DNA Supercoiling: An Ancestral Regulator of Gene Expression in Pathogenic Bacteria?[END_REF].

DNA supercoiling furthermore plays an important role in bacterial cells as an ancestral regulator of gene activity [START_REF] Dorman | DNA supercoiling is a fundamental regulatory principle in the control of bacterial gene expression[END_REF]. Indeed, as shown in figure 1A and 1B, high negative supercoiling levels, or hypercoiling (σ < σ basal ), favor higher expression of bacterial genes, as the thermodynamic reaction of promoter opening required to begin transcription is facilitated [START_REF] El Houdaigui | Bacterial genome architecture shapes global transcriptional regulation by DNA supercoiling[END_REF]. Conversely, high supercoiling levels, or DNA relaxation (σ > σ basal ), reduce gene expression. DNA hypercoiling and relaxation have been shown to affect the expression of 10% and 13% respectively of the genome of Streptococcus pneumoniae, in which supercoiling-sensitive genes can see up to 5-fold increases or decreases in their expression level when subjected to DNA relaxation (de la [START_REF] De La Campa | The Transcriptome of Streptococcus pneumoniae Induced by Local and Global Changes in Supercoiling[END_REF]. In Escherichia coli, DNA relaxation has been shown to affect the expression of around 7% of the genes [START_REF] Peter | Genomic transcriptional response to loss of chromosomal supercoiling in Escherichia coli[END_REF]. Moreover, in some bacterial species such as Buchnera aphidicola, an obligate aphid endosymbiont with a highly reduced genome [START_REF] Viñuelas | Conservation of the links between gene transcription and chromosomal organization in the highly reduced genome of Buchnera aphidicola[END_REF], gene regulation still takes place even in the near-total absence of traditional transcription factors; in these species, DNA supercoiling is suspected to be the main, if not the sole, regulatory mechanism [START_REF] Brinza | Genomic analysis of the regulatory elements and links with intrinsic DNA structural properties in the shrunken genome of Buchnera[END_REF]. When DNA is underwound (σ < σ basal , left), gene transcription rates are higher than when DNA is overwound (σ > σ basal , right). B. Promoter activity (equivalently, transcription level) e increases with the level of negative supercoiling -σ. C. The transcription of a gene by RNA polymerase (RNAP) generates a decrease in supercoiling upstream of the transcribed gene, and an increase downstream of the transcribed gene. D. Transcription-supercoiling coupling: the sign of the interaction between neighboring genes depends on their relative orientation.

Dynamic Properties of DNA Supercoiling

DNA supercoiling is under the influence of both internal and external constraints. It varies both in time, during the lifecycle of the bacterium, which alternates between growth and stationary phases [START_REF] Krogh | Impact of Chromosomal Architecture on the Function and Evolution of Bacterial Genomes[END_REF], and in space, as different regions of the chromosome experience different supercoiling levels [START_REF] Lal | Genome scale patterns of supercoiling in a bacterial chromosome[END_REF][START_REF] Junier | Conserved Units of Co-Expression in Bacterial Genomes: An Evolutionary Insight into Transcriptional Regulation[END_REF]. In bacteria, the supercoiling homeostasis is mainly regulated by the topoisomerases, but nucleoid-associated proteins (NAPs) such as H-NS also play a topological role by preventing the relaxation of superhelical stress at their fixation points in the genome, resulting in topological domains that have different supercoiling levels [START_REF] Krogh | Impact of Chromosomal Architecture on the Function and Evolution of Bacterial Genomes[END_REF].

Moreover, the global supercoiling level can change as a genome-wide response to extracellular stresses, such as changes in pH, osmolarity, temperature, or oxidative stress [START_REF] Duprey | The regulation of DNA supercoiling across evolution[END_REF].

Crucially, the transcription process itself plays a role in local supercoiling level variations. Indeed, when an RNA polymerase (RNAP in figure 1C) transcribes a gene, it follows the helical twist of the DNA template, but its rotation is hampered by frictional drag [START_REF] Ma | DNA supercoiling during transcription[END_REF]. The RNA polymerase therefore acts as a topological barrier and generates an accumulation of positive supercoiling through overwinding of the DNA molecule downstream of the transcribed region, and of negative supercoiling through underwinding upstream of the transcribed region, in what has been called the twin-domain model of DNA supercoiling [START_REF] Liu | Supercoiling of the DNA template during transcription[END_REF], as shown in figure 1C.

The Transcription-Supercoiling Coupling

As transcription is itself regulated by supercoiling (figures 1A and B), these interactions result in a possible dynamic coupling between the transcription levels of neighboring genes, which has been named the transcription-supercoiling coupling [START_REF] Martis | DNA Supercoiling: An Ancestral Regulator of Gene Expression in Pathogenic Bacteria?[END_REF]. As previously described, when a gene is transcribed (figure 1C), DNA upstream of the transcribed gene is underwound and the level of negative supercoiling increases (σ decreases), while DNA downstream of the transcribed gene is overwound and the level of negative supercoiling decreases (σ increases). This, in turn, affects the transcription of neighboring genes, as a higher level of negative supercoiling at the promoter of a gene increases its transcriptional activity. Therefore, a highly transcribed gene can increase the transcription level of the genes that are upstream of itself, and decrease the transcription level of the genes that are downstream of itself. The influence of this coupling on gene transcription is represented in figure 1D: when two neighboring genes lie in diverging orientations (figure 1D, left), the transcription of each gene generates a local increase in negative supercoiling around the other gene, thereby increasing the transcription level of that gene; each gene therefore reinforces the transcription of the other gene. Conversely, when two neighboring genes face each other in converging orientations (figure 1D, center), each gene is located downstream of the RNA polymerase during the transcription of the other gene, leading to a decrease in negative supercoiling and therefore a lower transcription level. In that case, each gene inhibits the other gene. Finally, if two genes are in a colinear orientation (figure 1D, right), the downstream gene up-regulates the upstream gene, and the upstream gene down-regulates its downstream neighbor.

Several mathematical and computational models have been proposed to describe the effect of the transcription-supercoiling coupling on the expression level of neighboring genes. In [START_REF] Meyer | Torsion-Mediated Interaction between Adjacent Genes[END_REF], a quantitative model of the supercoiling level at a locus of interest is proposed. DNA transcription is regulated by the opening energy of DNA around gene promoters, which directly depends on the supercoiling level. In this model, the reciprocal influence of neighboring genes can be obtained by computing the difference in transcription levels due to supercoiling and the subsequent variation in supercoiling, and iterating this system until a fixed point is reached. El [START_REF] El Houdaigui | Bacterial genome architecture shapes global transcriptional regulation by DNA supercoiling[END_REF] describe a more detailed stochastic model of DNA transcription involving explicit RNA polymerases and topoisomerases. The transcription level of a genomic region of interest is simulated using discrete time steps, during which RNA polymerases attach to the DNA template, progress along the transcribed region while generating positive supercoiling downstream and negative supercoiling upstream, and detach from the DNA, relaxing supercoiling constraints.

These models however limit themselves to mechanistic descriptions of the local interaction between genes, but do not try to generalize to the whole-genome scale nor to an evolutionary time frame. Yet, the dense gene content of bacteria suggests that the transcription-supercoiling coupling can generate a global transcriptional interaction network through the propagation of local supercoiling variations. Indeed, in bacteria, the distances between the beginning of two consecutive genes average around 1,000 bp [START_REF] Blattner | The Complete Genome Sequence of Escherichia coli K-12[END_REF]. This is low enough to connect multiple genes through the transcriptionsupercoiling coupling, as the typical distance at which this interaction operates is around a few thousand base pairs on each side of the transcribed gene [START_REF] El Hanafi | Activation and silencing of leu-500 promoter by transcription-induced DNA supercoiling in the Salmonella chromosome: Transcription-dependent modulation of leu-500 promoter in topA mutants[END_REF]. Measurements of the expression level of neighboring genes in bacteria have moreover experimentally demonstrated the existence of a coupling between transcription and supercoiling in specific gene systems, such as the ilvIH-leuO-leuABCD region in Escherichia coli or Salmonella typhimurium [START_REF] El Hanafi | Activation and silencing of leu-500 promoter by transcription-induced DNA supercoiling in the Salmonella chromosome: Transcription-dependent modulation of leu-500 promoter in topA mutants[END_REF][START_REF] Sobetzko | Transcription-coupled DNA supercoiling dictates the chromosomal arrangement of bacterial genes[END_REF][START_REF] Dorman | DNA supercoiling is a fundamental regulatory principle in the control of bacterial gene expression[END_REF]. Finally, experimental evidence of this coupling has also been obtained in a synthetic system using the ilvY and ilvC promoters of E. coli placed in diverging orientations on a plasmid, in which a decrease (resp. increase) in ilvY expression was correlated with a decrease (resp. increase) in ilvC expression [START_REF] Rhee | Transcriptional coupling between the divergent promoters of a prototypic LysR-type regulatory system, the ilvYC operon of Escherichia coli[END_REF].

In addition to a localized influence, global gene regulation through changes in the DNA supercoiling level has been shown to exist in nature. An example of this is Buchnera aphidicola, an endosymbiotic bacteria with a streamlined genome, in which control of the supercoiling level has evolved to be one of the main regulatory mechanisms [START_REF] Brinza | Genomic analysis of the regulatory elements and links with intrinsic DNA structural properties in the shrunken genome of Buchnera[END_REF]. Moreover, in Dickeya dadantii, a plant pathogenic bacteria, different genomic regions exhibit markedly different responses to changes in supercoiling [START_REF] Muskhelishvili | Coherent Domains of Transcription Coordinate Gene Expression During Bacterial Growth and Adaptation[END_REF], allowing the expression of pathogenic genes only in stressful environments. Finally, in the setting of experimental evolution, mutations in the regulation of supercoiling have been shown to drive the evolutionary response of Escherichia coli strains. In the Long-Term Evolution Experiment (LTEE) [START_REF] Lenski | Long-Term Experimental Evolution in Escherichia coli. I. Adaptation and Divergence During 2,000 Generations[END_REF], 12 populations of Escherichia coli have been maintained for over 80,000 generations, evolving and adapting to a glucose-limited environment.

In this experiment, parallel increases in the level of negative supercoiling been measured in 10 of the populations [START_REF] Crozat | Parallel Genetic and Phenotypic Evolution of DNA Superhelicity in Experimental Populations of Escherichia coli[END_REF], suggesting that a more negatively supercoiled genome confers an evolutionary advantage for that environment. Furthermore, mutations in two genes regulating DNA supercoiling, topA and fis, have been identified as the genetic basis for this phenotypic change, and have been verified to confer a fitness advantage [START_REF] Crozat | Long-Term Experimental Evolution in Escherichia coli . XII. DNA Topology as a Key Target of Selection[END_REF] when inserted into the genetic background of the ancestral strain. These results suggest a strong selection pressure to tune the level of DNA supercoiling to the new environment of the LTEE.

In conclusion, DNA supercoiling plays an important evolutionary role. It can serve as the support for the evolution of particular chromosomal organizations as a way to trigger certain sets of genes based on the change in supercoiling caused by specific environments, and its regulation is itself subject to selective pressure when adapting to new environments. The outer ring color shows locally high (σ > σ 0 , red) or low (σ < σ 0 , blue) supercoiling levels due to gene transcription. In the inner ring, closer genes interact more strongly (black arrows) than genes that are farther apart (gray arrows), either positively or negatively depending on their relative orientations.

Our model is written in Python and consists in an individual-based simulation, the source code of which is available at https://gitlab.inria.fr/tgrohens/evotsc/-/tree/alife-journal. It is also preserved for long-term archival [START_REF] Grohens | Archive of the EvoTSC repository on the Software Heritage archive[END_REF] using the Software Heritage archive [START_REF] Di Cosmo | Archiving and Referencing Source Code with Software Heritage[END_REF]). An individual in the model is represented by a circular genome (representative of most bacterial genomes), comprising a fixed number of genes, separated by non-coding intergenic regions.

Each gene is described by the following characteristics: its locus on the genome, its orientation, and its basal transcription (or expression) level. As we are mainly interested in the interplay between supercoiling and transcription, we voluntarily do not make the difference between gene expression levels, understood as mRNA or protein concentrations, and transcription levels, the immediate rate of mRNA production. Indeed, assuming a separation of timescales between the fast equilibrium of the transcription-supercoiling coupling, and the slow degradation of mRNAs, the concentration of a given mRNA is directly proportional to the transcription rate of its source gene.

Figure 2 shows the role played by the transcription-supercoiling coupling in an example genome. It includes the local supercoiling variations due to gene transcription, and the resulting gene interaction network, with each gene possibly activating or inhibiting its neighbors, depending on their relative orientations. Importantly to our approach, here genes do not interact only with their closest neighbors, but also with more distant genes, as is assumed to be the case in the gene-rich bacterial genomes (remember that E. coli genes are around one thousand base pairs apart, and that transcriptiongenerated supercoiling propagates around a few thousand base pairs on each side of the transcription site).

Mathematical Description of the Model

We model the transcription-supercoiling coupling between an individual's genes as a system of equations, which relate the supercoiling level at the locus of each gene σ i (for i ranging from 1 to n, the number of genes of the individual), and the expression level of every gene e i . The parameters of the system are described by the genome of the individual, as will be detailed below.

In our model, the supercoiling at a given locus on the genome depends on three factors: the individual's basal supercoiling level σ basal , the variation in supercoiling due to environmental conditions σ env , and the variation in supercoiling due to the transcription of the neighboring genes. We compute this local variation in supercoiling at the locus of each gene with the help of a gene interaction matrix, whose coefficient at position (i, j) describes the influence of gene j on gene i. The coefficients are given by the following equation:

∂σ i ∂e j = η • c • max(1 - d(i, j) d max , 0) (1)
More precisely, the interaction level between two genes depends on the relative orientation of the genes, as the transcription of a gene increases supercoiling at the locus of downstream genes and decreases supercoiling at the locus of upstream genes (remember that an increase in supercoiling means a decrease in transcription). Therefore, we choose η = 1 if gene i is downstream of gene j and η = -1 otherwise (if i = j, η = 0 as a gene does not interact with itself). The interaction level also depends on gene distance, as genes that are further apart on the genome interact less strongly, so the strength of the interaction linearly decreases with the intergenic distance d(i, j), and reaches 0 when d(i, j) = d max , the maximum distance above which the interaction vanishes. Finally, an interaction coefficient c is applied to adjust the strength of the coupling.

Using this interaction matrix, we compute the level of supercoiling σ i at the locus of every gene, which depends on the transcription level of all the other genes, on the basal supercoiling and on the environmental supercoiling:

σ i = σ basal + σ env + n j=1 ∂σ i ∂e j e j (2)
The transcription level e i of every gene as a function of total supercoiling is then modeled with a sigmoidal activation curve, following El Houdaigui et al. ( 2019). The equation is given by:

e i = 1 1 + e (σi-σ0)/ε
(3)

In this equation, σ 0 is a parameter that represents the inflexion point of the sigmoid, that is the supercoiling level at which the gene is at half its maximum transcription rate, and ε a scaling factor that represents the strength of the dependence of the transcription level on the supercoiling level.

Finally, in order to obtain the phenotype of an individual, we numerically compute a solution to the system of equations 2 and 3, using a fixed point algorithm. This solution represents the state (of gene expression and supercoiling at every locus) towards which the individual would converge over time. Let f (e i ) be the function that computes new supercoiling levels σ ′ i from e i using equation 2, then computes new expression levels e ′ i from the new σ ′ i using equation 3, and finally returns e ′ i . In order to compute a fixed point of f , that is a set of transcription levels e * i such that f (e * i ) = e * i , we start with the basal transcription levels e 0 i that are a property of each gene, and iterate the sequence

e t+1 i = 1 2 (e t i + f (e t i )
), until the difference between two successive iterations is below a given threshold.

In our setting, this algorithm has empirically always converged to a solution that is a stable fixed point of the function, and that is therefore interpretable from a biological perspective.

10 Figure 3 shows the genome (left, outer ring) of an example individual with a genome of 13,000 bp and n = 13 genes evenly spaced along the genome, and with a basal supercoiling of σ basal = -0.06.

The basal transcription level of each gene is randomly chosen between 0 and 1, and the iterations of the fixed point algorithm giving the final gene transcription levels are shown on the right. In this individual, the non-linear effect of the interaction between neighboring genes is clearly visible. Indeed, six genes (A, B, D, E, H, and I) end up at a high transcription rate at the fixed point (or solution) of the system, while the others end up at low transcription rates. These activated genes can be grouped into 3 pairs (A and B, D and E, H and I), all of which are pairs of adjacent genes in divergent orientations. Even though gene D has a low (around 0.3) basal transcription rate, it eventually reaches a high transcription state because of its positive interaction with gene E. Conversely, genes F and G start with a high transcription rate, but are repressed by their neighbors H and E, and are therefore silenced when the system converges. We can also observe complex behaviors in the model, as the gene expression levels pass through very different states during convergence to the solution. Indeed, the transcription level of gene K initially increases due to its interaction with gene J, but both genes end up in a low transcription state, as they are inhibited by the very active gene I. The final supercoiling level along the genome (left, inner ring) moreover demonstrates the effect of the transcription-supercoiling coupling on local supercoiling. Highly transcribed genes, such as A and B, generate a large variation in the supercoiling level on their upstream and downstream sides, and the positive feedback loop between genes in divergent pairs is made clear by the very high negative value of the supercoiling level between each gene in these two pairs. DNA is severely overwound compared to normal, no gene at all is activated (with an expression level e > 0.5). As the external influence of the environment on supercoiling decreases to σ env = 0, corresponding to normal relaxation of DNA, and then to σ env = -0.1, 8 of the 13 genes of the individual reach an activated state. Finally, for σ env = -0.2, there is a strong environmental pressure towards high gene transcription levels, and most genes are indeed activated; however, even at this level of σ env , some genes remain shut down, because of the high amount of positive supercoiling (in red) generated by the transcription of their neighbors. reversing the positions and orientations of genes B and C. This is an example of a genomic inversion, which will be presented in further detail in section 3.3. The start point of this inversion falls between genes A and B and its end point between genes C and D; this results in the reversal of segment [BC] relative to the rest of the genome. Here, we can see that the diverging orientation that was present between genes A and B has vanished, replaced by a set of genes in colinear orientation, from A to D.

Effect of the Environmental Supercoiling on Gene Activation Levels

Influence of Relative Gene Positions on Gene Activation Levels

This genomic reorganization results in the loss of the activation of genes A and B, as gene B is now more strongly inhibited by gene D due to its closer genomic location, and as genes A and B are not in a positive feedback loop due to diverging orientations any longer; only the pairs of genes D and E, and H and I, remain activated.

Based on these observations, we can confirm that in our model, the transcription-supercoiling coupling generates complex networks of genome-wide interactions between genes, and that these networks are directly dependent on the architecture of the genome.

3 An Evolutionary Genome-Wide Model of the Transcription-

Supercoiling Coupling

Having shown that transcriptional activity depends on the organization of the genome, we now question to which extent evolution can simultaneously leverage the organization of the genome and the transcription-supercoiling coupling in order to adapt gene regulatory activity to different environments.

Indeed, as has been observed in Dickeya dadantii [START_REF] Muskhelishvili | Coherent Domains of Transcription Coordinate Gene Expression During Bacterial Growth and Adaptation[END_REF], different phenotypes can evolve as a response to different supercoiling levels induced by the environment, and the transcriptionsupercoiling coupling could play a role in enabling the existence of this reaction norm.

In this section, we expand our model into an evolutionary simulation. At each generation of the simulation, all individuals are evaluated and their fitness values are computed, based on their gene transcription levels. Then, the individuals of the new generation are chosen by picking their ancestor from the current generation, with a probability proportional to the ancestor's fitness. The model is panmictic, meaning that any individual in the population can be chosen as the ancestor of any new individual. Finally, during replication, the genome of each new individual stochastically undergoes a number of mutations, before the new individual is evaluated again; importantly, these mutations do not impact genes themselves, but only the spatial organization of the genome: gene orientations, syntenies, and intergenic distances.

Evolutionary Model: Evolution in Two Separate Environments

We model the evolution of populations of individuals that experience two different environments, named A and B. Each environment is defined by its value of σ env , respectively σ A and σ B , which represent the change in the supercoiling level due to the environment [START_REF] Dorman | DNA supercoiling is a fundamental regulatory principle in the control of bacterial gene expression[END_REF].

In order to have environments with distinct effects, we choose a value of σ A = 0.1, for which isolated genes are effectively inhibited (as in the top-left panel of figure 4), and a value of σ B = -0.1, for which some but not all genes are activated (bottom-left panel).

We separate genes into three classes, based on the environments in which they must be activated:

either in both environment A and environment B (AB genes), only in environment A (A genes), or only in environment B (B genes). These classes allow us to define optimal phenotypes for both environments: in environment A, both A and AB genes should be activated, whereas B genes should be inhibited. Conversely, in environment B, only B and AB genes should be activated, but not A genes.

Fitness

In order to compute the fitness of an individual, we define an optimal phenotype ẽA (resp. ẽB ), corresponding to the vector of the expected expression level ẽA i for each gene i in environment A (resp. environment B). We choose an expected expression level of ẽ = 1 for genes that should be activated, which corresponds to the maximum possible expression level of a gene in our model. Similarly, we choose ẽ = 0 for genes that should be inhibited, which is the minimum expression level that is attainable. Then, in each environment, we compute the gap g A (resp. g B ), or average square distance of the individual's gene transcription levels e A (the vector constituted by the transcription level e A i of each gene i) to the optimal levels ẽA (resp. e B and ẽB ). The gap g A is computed as follows:

g A (e A ) = 1 n n i=1 (e A i -ẽA i ) 2 (4)
The gap g B is computed in the same way. Finally, we compute the fitness of the individual by summing the gap in each environment, and applying an exponential scaling: f = e -k(g A +g B ) , where k is a scaling factor representing the selection pressure. A higher value of k means that well-adapted individuals, those which have a smaller gap, will have an even higher fitness value compared to other individuals; we typically use k = 50, meaning that a small decrease in the gap compared to other individuals yields a large reproductive advantage.

Mutational Operator: Genomic Inversions

We introduce only one kind of mutation in our model, which is genomic inversions: we choose two breakpoints randomly on the genome, and reverse the genomic content between these points. Genes are then reinserted in the genome in the opposite orientation and order, taking care to update all intergenic distances appropriately. Note that in our model, genes have a length of zero and the break points can therefore not fall inside a gene. Moreover, an inversion has no effect if both breakpoints fall between two neighboring genes (as only an intergenic region would be affected), but can impact any number of genes otherwise. Genomic inversions hence affect gene syntenies and orientations, and As a result of this inversion, genes B and C have been switched from the forward to the backward orientation, and the intergenic distances between A and C on the one hand, and B and D on the other hand, have been modified; however, the relative orientation of B and C, and hence their interaction subnetwork, remain unchanged. This results in changes to the gene interaction network: instead of mutual activation between genes A and B and mutual inhibition between genes C and D, all four of those genes now lie in colinear orientations, in which each of these genes activates its upstream neighbor but represses its downstream neighbor.

Experimental Setup and Parameter Values

We initialized the simulation with a clonal population of N = 100 copies of an initial individual with the following genome: 60 genes in random orientations, uniformly distributed along a 60,000 bp genome, and equally divided between the AB, A and B classes. We chose a maximum interaction distance of d max = 2500, meaning that each gene initially interacts with its 2 closest neighbors in each direction through the transcription-supercoiling coupling. Note that as inversions may change intergenic distances, genes can become closer or further apart during evolution. We set the basal supercoiling level σ basal to the average supercoiling level in E. coli of -0.06 [START_REF] Crozat | Long-Term Experimental Evolution in Escherichia coli . XII. DNA Topology as a Key Target of Selection[END_REF], and σ 0 to -0.06 as well, so that in the absence of other sources of supercoiling (either environmental or through the coupling), the default activity level of a gene is 0.5. Finally, we set c = 0.3, in order to have comparable values for the variations in supercoiling due to the environment and due to the transcription-supercoiling coupling, and ε = 0.03, so that the variations in supercoiling have a qualitatively mild effect on gene expression.

In order to run the simulations, we evolved 15 different populations for 250,000 generations; the simulation lasted for approximately 48h on a computer with Intel Xeon E5-2640 v3 @ 2.60GHz CPUs, using around 100 MB of RAM per replicate.

Adaptation of Gene Expression Levels to Different Environments

Figure 7 summarizes the differences in the proportion of activated genes for each of the three sets of genes, between environments A and B, averaged over the 15 repetitions. In the figure, we consider a gene to be activated if its activity at the end of the lifecycle is over 0.5, and look at the average proportion of activated genes in the best individual of every replica. Recall that the evolutionary target for AB genes is an expression level of 1 in both environments, for A genes an expression level of 1 in environment A and 0 in B, and vice-versa for B genes. After 250,000 generations of evolution, individuals have acquired genomes that allow all AB genes to be activated in both environments, and that allow all B genes to be activated in environment B and inhibited in environment A. On average, over 60% of A genes are activated in environment A, which imposes a positive change in supercoiling (σ A = 0.1) and makes gene activation harder. Conversely, less than 5% of A genes are activated in environment B, in which it is easier for genes to be activated (σ B = -0.1). The final expression levels of A genes therefore show that specific sets of genes can be activated by the transcription-supercoiling coupling despite environmental hurdles. Finally, details of the evolution of one of the 15 replicate populations are shown in figure 9. We can first see that the number of activated AB genes of the best individual at each generation quickly rises to 20 (out of 20 genes of that type) in both environment A and environment B; this shows that evolving a phenotype that is resistant to environmental perturbations, having genes that are always activated, is easy in the model. For A genes and B genes, we observe an asymmetric tendency during the course of evolution towards activation in the target environment, and inhibition in the opposite environment. However, the difference in the number of activated B genes between environment A and environment B is much higher than for A genes. As already mentioned above, this asymmetry comes from the different requirements expected of A genes and B genes: gene activation is easier in environment B than in environment A, as it is easier for a gene to become activated in an environment with a lower overall supercoiling level. A genes therefore have to be activated in a harder environment, and inhibited in a simpler environment, whereas B genes have to do the opposite.

Discussion and Perspectives

DNA supercoiling plays a fundamental role in the regulation of gene transcription in bacteria, and an important part of this role could be mediated by the local variations in supercoiling caused by the transcription-supercoiling coupling. While the influence of the global supercoiling level on gene transcription [START_REF] Lal | Genome scale patterns of supercoiling in a bacterial chromosome[END_REF][START_REF] Ma | DNA supercoiling during transcription[END_REF][START_REF] Dorman | DNA supercoiling is a fundamental regulatory principle in the control of bacterial gene expression[END_REF][START_REF] Martis | DNA Supercoiling: An Ancestral Regulator of Gene Expression in Pathogenic Bacteria?[END_REF], the evolutionary importance of supercoiling regulation [START_REF] Crozat | Long-Term Experimental Evolution in Escherichia coli . XII. DNA Topology as a Key Target of Selection[END_REF][START_REF] Crozat | Parallel Genetic and Phenotypic Evolution of DNA Superhelicity in Experimental Populations of Escherichia coli[END_REF][START_REF] Duprey | The regulation of DNA supercoiling across evolution[END_REF] and the mechanistic details of the transcription-supercoiling coupling [START_REF] Meyer | Torsion-Mediated Interaction between Adjacent Genes[END_REF][START_REF] El Houdaigui | Bacterial genome architecture shapes global transcriptional regulation by DNA supercoiling[END_REF] have all already been studied, no existing work did to our knowledge tackle the question of the possible role of the transcription-supercoiling coupling at both the whole-genome scale and on an evolutionary time scale.

In this work, we have developed a genome-wide model of the influence of DNA supercoiling on gene transcription, incorporating both the global influence of the environment and the local variations due to the transcription-supercoiling coupling on the supercoiling level. We have shown that, in our model, complex interactions implicating several genes emerge from the coupling between supercoiling and transcription. Indeed, A genes display an activation pattern that would not be obtainable without the network of interactions that results from the coupling. Thanks to this network, A genes are activated in an environment where isolated genes would be inhibited, and inhibited in an environment where isolated genes would be activated. The transcription-supercoiling coupling therefore enables the selective activation or inhibition of specific sets of genes, providing a non-monotonic response to environmental variations through changes in the level of DNA supercoiling. Furthermore, we have shown, using an in silico experimental evolution approach, that natural selection can leverage this biophysical mechanism to selectively turn on or off several pools of genes, using only the very simple mutation operator of genomic inversions, that affect the relative positions and orientations of genes on the genome but do not change genome length or basal gene transcription rates, and that this behavior is able to evolve under a wide range of parameter values. This response of gene transcription levels to DNA supercoiling reflects a phenomenon which has been observed in vivo in the expression of pathogenicity-related genes in specific environments, such as the normally lethal inside of the macrophage for the mammalian pathogen S. enterica [START_REF] Cameron | Transmission of an Oxygen Availability Signal at the Salmonella enterica Serovar Typhimurium fis Promoter[END_REF], or plant tissue for D.

dadantii [START_REF] Hérault | Role of the LysR-type transcriptional regulator PecT and DNA supercoiling in the thermoregulation of pel genes, the major virulence factors in Dickeya dadantii: Dickeya dadantii PecT protein and virulence thermoregulation[END_REF].

Our model voluntarily stays very simple, only incorporating the most important feature of the transcription-supercoiling coupling, the non-linear interaction between the expression levels of neigh-boring genes. This simplicity therefore hints at the possible pervasiveness of this regulation mechanism throughout the prokaryotic realm. Nonetheless, in order to go further and represent more accurately the diversity of gene behaviors found in real life, several more dimensions could be integrated to the model. At present, the target for genes in our model is bistability, meaning that genes should end up fully activated or fully inhibited. A more biologically plausible approach would be to relax this restriction and give genes arbitrary expression targets, in order to determine to which extent the transcription-supercoiling coupling is able to finely regulate gene expression. Furthermore, unlike in our model (in which all genes have the same response curve to DNA supercoiling), the genes of biological organisms can show different responses to the supercoiling level. These differences are in part caused by the GC content at the gene promoter [START_REF] Forquet | Role of the Discriminator Sequence in the Supercoiling Sensitivity of Bacterial Promoters[END_REF], and some genes can even respond in the opposite direction to changes in the supercoiling level, that is to say be activated rather than inhibited by DNA relaxation. This behavior is for instance present in the gyrA and gyrB genes that encode the gyrase subunits in E. coli [START_REF] Peter | Genomic transcriptional response to loss of chromosomal supercoiling in Escherichia coli[END_REF]. Moreover, while our model studies its role in an abstract transcription model, supercoiling intervenes during different parts of transcription initiation, and in transcript elongation and transcription termination as well [START_REF] Martis | DNA Supercoiling: An Ancestral Regulator of Gene Expression in Pathogenic Bacteria?[END_REF].

Incorporating such precise mechanistic processes into our model could give more accurate information on the link between the position of genes on the genome and their transcription rate. Similarly, increasing the number of genes of individuals in our model to match bacterial gene numbers might provide more fine-grained results, but is computationally intractable in the current implementation of the model. Furthermore, investigating the behaviors of individuals when they are placed successively in different environments, rather than evaluated separately in each environment, would also bring more information on the plasticity of the network of gene interaction levels that emerges from the transcription-supercoiling coupling. Finally, another valuable approach in order to bring this model closer to biology would be to incorporate it into a larger existing framework, such as the Aevol in silico experimental evolution platform [START_REF] Rutten | Adapting the engine to the fuel: Mutator populations can reduce the mutational load by reorganizing their genome structure[END_REF], which models the bacterial genome in much more detail, in order to leverage the power of a well-understood digital organism model.

  Figure1: A. When DNA is underwound (σ < σ basal , left), gene transcription rates are higher than when DNA is overwound (σ > σ basal , right). B. Promoter activity (equivalently, transcription level) e increases with the level of negative supercoiling -σ. C. The transcription of a gene by RNA polymerase (RNAP) generates a decrease in supercoiling upstream of the transcribed gene, and an increase downstream of the transcribed gene. D. Transcription-supercoiling coupling: the sign of the interaction between neighboring genes depends on their relative orientation.

  Both the importance of supercoiling regulation and the detailed mechanisms of the transcriptionsupercoiling coupling at the local scale have well been studied, but a thorough analysis of the genomewide effect of the transcription-supercoiling coupling on gene expression, and of its possible evolutionary use by natural selection, remains missing. In this work, we describe a new model which incorporates a high-level model of global supercoiling regulation and of the transcription-supercoiling coupling within an in silico experimental evolution setting. Using this model, we first investigate the non-linear variation in gene transcription levels at the whole-genome scale in response to variations in the global supercoiling level. Then, we study the evolutionary trajectory of gene activation patterns in individuals subjected to different environments.We show that in our model, a genome-scale gene interaction network emerges from local interactions, and creates a reaction norm in response to the change of a single parameter, the global supercoiling level, caused by different environments. Moreover, we demonstrate that, using genomic inversions as the only mutation operator, and therefore only changing the relative positions and orientations of genes on the genome, evolution can select genomes displaying qualitatively different phenotypes in different environments, characterized by different global supercoiling conditions. This article is an extended version of[START_REF] Grohens | A Genome-Wide Evolutionary Simulation of the Transcription-Supercoiling Coupling[END_REF].

Figure 2 :

 2 Figure 2: Genes along an example genome and local variations in supercoiling (outer ring), and the associated gene interaction network (inner ring). The outer ring color shows locally high (σ > σ 0 , red) or low (σ < σ 0 , blue) supercoiling levels due to gene transcription. In the inner ring, closer genes interact more strongly (black arrows) than genes that are farther apart (gray arrows), either positively or negatively depending on their relative orientations.

Figure 3 :

 3 Figure 3: Left: genome (outer ring) and stable state level of supercoiling σ (inner ring) of an example individual with 13 genes in the model. Right: transcription levels of the individual's genes during the iterations of the fixed point computation, in an environment given by σ env = 0.05. Solid lines represent genes in forward orientation, and dashed lines genes in reverse orientation.

Figure 4 :

 4 Figure4: Influence of the environment supercoiling σ env on the stable state local supercoiling level (top row) and gene transcription levels (bottom rows) of the example individual. From left to right and top to bottom: at σ env = 0.1, no genes are activated (e > 0.5); at σ env = 0.0 and at σ env = -0.1, 8 genes are activated; at σ env = -0.2, 10 genes are activated. Lower values of σ env result in the activation of more genes, reflecting the in vivo effect of higher negative supercoiling.

Figure 4

 4 Figure 4 captures the influence of the environmental change in supercoiling σ env on the local supercoiling level due to the transcription-supercoiling coupling (top row) and on the repartition of genes between the activated and inhibited states (bottom rows), again using the example individual already shown in figure 3. From left to right and top to bottom: at a high value of σ env = 0.1, meaning that

Figure 5 :

 5 Figure 5: Genome, local supercoiling and gene expression levels of a new individual obtained from the individual in figure 3 by switching the positions and orientations of genes B and C.

Figure 5

 5 Figure 5 again shows the local supercoiling and gene expression levels of the individual in figure 3, after

Figure 6 :

 6 Figure 6: Result of the inversion of a genomic segment containing genes B and C from the individual presented in figure 2. The gene interactions which have changed due to the inversion are drawn in red. This illustration genome corresponds to the actual individual in our model presented in figure 5.

Figure 6

 6 Figure 6 presents a genome obtained by performing an inversion on the genome shown in figure 2.

Figure 7 :FitnessFigure 8 :

 78 Figure7: Fraction of activated genes of each type in each environment at the end of the lifecycle, averaged over the best individual in the last generation of each replica. Boxplots represent the median and quartiles, and dots flier data points. For A genes and B genes, activation levels differ depending on the environment: p-value 2.40 × 10 -17 for A genes, and p-value < 1 × 10 -25 for B genes (Student's t-test for dependent samples).

Figure 9 :Figure 10 :

 910 Figure 9: Number of activated genes of each type and fitness of the best individual at every generation of replicate 13, with a population size of N = 100, for 250,000 generations. The number of active AB genes increases until it reaches 20, in both environment A (top) and environment B (bottom). The number of active A (resp. B ) genes increases in environment A (resp. B) and decreases in environment B (resp. A) over time, thus converging towards their evolutionary target.

Our experimental results show that, in a model of gene transcription that is structured around the transcription-supercoiling coupling, complex gene interaction networks can evolve. These gene interaction networks are sensitive to environmental variations, which are mediated in our model by a single parameter: σ env , the amount of global supercoiling that is due to the environment.

Robustness of Gene Network Evolution

In order to ensure that our results remain experimentally valid over a broad range of parameter values, we ran additional sets of simulations, in which we changed respectively the sensitivity of gene promoters to supercoiling changes (ε in equation 3), the interaction coefficient used in computing the local supercoiling due to the transcription-supercoiling coupling (c in equation 1), and the strength of the change in supercoiling imposed by the environment (σ A and σ B ). We chose sets of logarithmically-spaced values for each parameter, and ran 5 replicates of the evolution experiment for 250,000 generations for each parameter value. Note that, for extreme parameter values, gene expression levels did in some cases not converge to stable states by the maximum number of computation steps. In this situation, we chose to retain the gene expression levels at the last step as the phenotype of the affected individuals. For all values of ε except 0.1, the behavior from the main run is qualitatively replicated.

The results of the additional simulations are presented in figures 11, 12 and 13. For ε, we chose values of ε = 0.003, ε = 0.01, and ε = 0.1, compared to an initial value of ε = 0.03, and the results are shown in figure 11. For the values of ε lower than the default (top row), representing a higher sensitivity of promoters to supercoiling, we observe the evolution of differentiated gene expression levels as in the main run (bottom-left panel), whereas for the higher value of ε (bottom-right panel), A genes are not expressed in environment A by the end of evolution. In this case, promoters are not sensitive enough to the supercoiling variations caused by the transcription-supercoiling coupling, and genes are unable to overcome the highly positive supercoiling of environment A. Finally, we also investigated different amplitudes in the difference in supercoiling level between the two environments, by choosing values of σ A = 0.01, σ A = 0.05 and σ A = 0.2, and σ B = -σ A respectively in each case (for an initial value of σ A = 0.1 and σ B = -0.1). We observe that, when σ A = 0.2 (bottom-right panel), the environmental supercoiling constraint is too high and A genes are not activated in environment A by the end of the runs. However, for environments closer to each other than the default (top row), evolution is able to leverage the differences in supercoiling between these environments to evolve differentiated phenotypes, as in the main run (bottom-left panel), showing that our model remains sensitive to small changes in environmental supercoiling.

To conclude, in our model, the gene interaction network is therefore able to respond to different environments and can evolve an efficient regulation of gene expression under a broad range of parameter values, reinforcing the hypothesis that a supercoiling-mediated coupling between gene expression levels could indeed play a functional role in biological organisms.


Antoine Scherrer

Antoine Fraboulet

Tanguy Risset

Long-Range Dependence and On-chip Processor Traffic

Keywords: Network on Chip, System on Chip, embedded software, long range dependence, network traffic

HAL is

Introduction

The next generation of embedded systems on chip (SoC) will soon integrate many processor cores together with dedicated circuits connected via a network. This is particularly true for chips integrated in telecommunication and multimedia devices which require more and more computing power. These devices must be designed faster to meet the market demand which implies that more and more embedded software is present on these chips. Embedded software brings a faster and more flexible design process. However, as performance is decreased compared to dedicated hardware, parallelization and efficient communication schemes must be implemented to reach the required computing power. In particular, controlling network throughput and contention is essential and is in general very difficult, because of non determinism introduced by parallelism. This is why network prototyping will become a major issue in future SoC design processes.

Recent advances in Internet traffic analysis have shown that predicting performance of a network is a difficult task. Stochastic modeling of the traffic must be used and it has been shown [START_REF] Paxon | Wide-area traffic: The failure of poisson modeling[END_REF] that this modeling must take into account second order stochastic properties (covariance in addition to marginal law). More recently Varatkar and Marculescu [START_REF] Varatkar | On-chip traffic modeling and synthesis for mpeg-2 video applications[END_REF] have shown that on-chip traffic might require a similar modeling. They have found long-range dependent behavior in communications between different parts of the Mpeg-2 video decoding application. The presence of long-range dependence (lrd) in embedded applications, as in Internet traffic, could imply that ips communication behavior will not be correctly modeled by a simple random traffic even if it is adjusted to a particular first order statistical law. However, the experiments of Varatkar and Marculescu were not done at cycle accurate level and it is still not clear how the lrd will impact the network-on-chip.

We want to emphasize the originality of our work with respect to the seminal paper of Varatkar and Marculescu [START_REF] Varatkar | On-chip traffic modeling and synthesis for mpeg-2 video applications[END_REF], in which it was shown that the aggregated throughput of the macro-block arrival process at the IDCT/IQ module of a hardware Mpeg-2 decoder exhibits long-range dependence. Standard queuing analysis with self-similar input was then used in order to evaluate the buffer usage. In the present work, we first investigate the presence of longrange dependence in the traffic produced by an on-chip processor ip at the cycle accurate level and then we show that lrd has a very different impact when the low level communication protocol uses non-split transactions. In practice, if the application is executed by a processor connected to a cache, the fact that each read transaction is waiting for an answer (this is called a non-split transaction), makes the impact of lrd negligible. If the application is executed by a dedicated ip that computes on data streams (i.e. not waiting for answers: split transactions), lrd decreases performance.

We also demonstrate here a way to adjust the size of the fifos given particular statistical properties of the executed application. As pointed out in [START_REF] Bjerregaard | A survey of research and practices of networkon-chip[END_REF], "buffers account for the main part of the NoC router area, it is a major concern to minimize the amount of buffering necessary under performance requirements". We quantify experimentally the impact of lrd on fifo usage and on average communication latency and hence propose a way to evaluate the required size of fifos to cope with the presence of lrd. This is possible because we developed a complete NoC prototyping environment [START_REF] Scherrer | Automatic phase detection for stochastic on-chip traffic generation[END_REF][START_REF] Scherrer | Generic multi-phase on-chip traffic generator[END_REF] in which we are able to produce lrd traffic. The NoC used in the simulations is dspin [START_REF] Greiner | A generic architecture for on-chip paquets-switched interconnections[END_REF], and the simulation is done at cycle accurate level: routers are precisely simulated and instrumented in order to observe the exact behavior of what would occur on a real chip.

The paper is organized as follows: the next section presents the theoretical background needed to understand the notion of lrd, Section 3 presents our NoC prototyping environment. Experiments and results about the presence and the impact of lrd are presented in Section 4. Related works are referenced in Section 5.

Network traffic and long-range dependence

We start this section by introducing a simple on-chip traffic modeling formalism. We then give theoretical background on long-range dependence whereof the impact on performance of on-chip networks is the main topic of the paper. The traffic produced by a component is modeled as a sequence of transactions composed of flits (flow transfer units) corresponding to a bus-word. The k th transaction is a 5-uple (A(k), C(k), S(k), D(k), I(k)) meaning respectively, target address, command (read or write), size of transaction, delay and interrequest time. This is illustrated in Fig. 1. We also define the latency of the k th transaction L(k) as the number of cycles between the start of the k th request and the start of the associated response. This is basically the roundtrip time in the network and will be used in our experiments to evaluate contention. One can distinguish two main communication schemes used by ips: the non-split transactions scheme where the ip is not able to send a request until the response to the previous one has been received, and the split transactions scheme in which new requests can be sent without waiting for the responses. The non-split transaction scheme is used by processors and caches for instance (although for cache, it might depend on the cache parameters). The split transaction scheme is used by dedicated ips performing computation of streams of data which are transmitted via dedicated direct memory access (dma) modules.

On-chip traffic modeling

L(k+1) R(k) S(k+1) R(k+1) A(k),C(k) A(k+1),C (k+1)
Our goal is to emulate real traffic by means of traffic generators. With the above formalism, traffic generation means producing (A(k), C(k), S(k), D(k), I(k)), for each k. This can be done either deterministically (replay of a recorded trace), or using stochastic processes. In this work, we use stochastic traffic generators. Note that, in this case, we make the assumption that the elements of the transaction sequence are independent.

Let us recall that a stochastic process X is a sequence of random variables X[i] (we use brackets to denote random variables). We will consider two statistical characteristics of stochastic processes: the marginal law, which represents how the values taken by the process are distributed, and the covariance function, which gives information on the correlation between the random variables of the process as a function of the time lag between them. For instance, the sequence of delays D(k) can be generated as the sample path of some stochastic process {D[i]} i∈N , with prescribed first and second statistical orders.

Long-range dependence

Definitions

Long-range dependence (lrd) is a property of a stochastic process which is defined as a slow decrease of its covariance function [START_REF]Self-Similar Network Traffic and Performance Evaluation[END_REF]. This function describes how the random variables of a process are correlated with each-other as a function of the time lag between these random variables. The covariance function γ X of a stochastic process {X[i]} i∈N is defined as follows (E is the expectation):

γ X (i, j) = E(X[i]X[j]) -E(X[i])E(X[j])
If the process X is wide-sense stationary, then its mean is constant (

∀(i, j) ∈ N 2 , E(X[i]) = E(X[j]) ∆ = E(X)
) and its covariance reduces to a one variable function as :

∀(i, j) ∈ N 2 , γ X (i, j) = γ X (0, |i -j|) ∆ = γ X (|i -j|)
We expect this function to be decreasing, because correlated data are more likely to be close (in time) to each other. However, if the process is long-range dependent, then the covariance decays very slowly and is not summable:

k∈N γ X (k) = ∞
Therefore, long-range dependence reflects the ability of the process to be highly correlated with its past, because even at large lags, the covariance function is not negligible. This property is also linked to self-similarity which is more general: in can be shown that asymptotic second order self-similarity implies long-range dependence [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF].

lrd is an ubiquitous property of Internet traffic [START_REF] Leland | On the self-similar nature of ethernet traffic (extended version)[END_REF][START_REF] Paxon | Wide-area traffic: The failure of poisson modeling[END_REF], and it has also been attested on on-chip multimedia (MPEG-2) applications by Varatkar and Marculescu [START_REF] Varatkar | On-chip traffic modeling and synthesis for mpeg-2 video applications[END_REF]. The main interest in lrd resides in its strong impact on network performance [START_REF]Self-Similar Network Traffic and Performance Evaluation[END_REF]. Above all, the needed memorization in the buffers is higher when the input traffic has this property [START_REF] Erramilli | Experimental queueing analysis with long-range dependent packet traffic[END_REF]. As a consequence, for macronetworks as well as for on-chip networks, lrd should be taken into account if it is manifested in the traffic that the network will have to handle.

A long-range dependent process is usually modeled with a power-law decrease of the covariance function as follows :

γ X (k) ∼ k→+∞ ck -α
The exponent α (also called scaling index ) provides a parameter to quantify the process' long-range dependency (0 < α ≤ 1). The Hurst exponent, denoted H, is the classical parameter for describing self-similarity [START_REF]Self-Similar Network Traffic and Performance Evaluation[END_REF]. Because of the analogy between lrd and self-similarity, it can be shown that a simple relation exists between H and α: H = (2 -α)/2. As a consequence, H (1/2 < H < 1) is the commonly used parameter for characterizing long-range dependence. Note that when H = 0.5, then there is no long-range dependence (this is also referred to as short-range dependence).

Estimation of the Hurst parameter

We have used a standard wavelet-based methodology for the estimation of the Hurst parameter [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF] in our experiments. Let ψ j,k (t) = 2 -j/2 ψ 0 (2 -j t -k) denote an orthonormal wavelet basis, designed from the mother wavelet ψ 0 . The j index represents the scale: the larger j is, the more the wavelet is dilated. The k index is a shift in time.

For any (j, k), d X (j, k) = ψ j,k , X are called the wavelet coefficients of the stochastic process X ( ., . is the inner product in the L 2 functional space).

These wavelet coefficients permit a study of the X process at various times (values of k) and various scales (values of j).

When X is a long-range dependent process with parameter H, the following limit behavior for the expectation of wavelet coeficients can be shown [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF]:

∀j, E(d X (j, k) 2 ) ∼ j→+∞ C2 j(2H-1) (1)
It can also be shown that the time averages S j for each scale j (n j is the number of wavelet coefficients available at scale j):

S j = (1/n j ) n j k=1 |d X (j, k)| 2 (2)
can be used as relevant, efficient and robust estimators for E(d X (j, k) 2 ) [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF]. From Eq. ( 1) and ( 2), the estimation of H is as follows: i) plot log 2 S j versus log 2 2 j = j and ii) perform a weighted linear regression in the limit of the coarsest scales (see for instance Fig. 2). These plots are commonly referred to as log-scale diagrams (LD) . In such diagrams, lrd is demonstrated by a straight line behavior in the limit of large scales. In particular, if the line is horizontal, then H = 0.5 and there is no long-range dependence.

To illustrate how we use this tool to evaluate the Hurst parameter, we provide in Fig. 2 a typical log-scale diagram extracted from an Internet trace [START_REF]Wand network research group[END_REF].

Along the x axis are the different values of the scale j at which the process is observed. For each scale, log 2 S j is plotted together with its confidence interval (vertical bars). The Hurst parameter can be estimated if the different points plotted are aligned on a straight line for large scales.

Synthesis of long-range dependent processes

The synthesis (generation of sample paths) of long-range dependent processes is easy if the marginal law is Gaussian [START_REF] Bardet | Long-Range Dependence: Theory and Applications, chap. Generators of long-range dependent processes: A survey[END_REF]. The so-called Fractional Gaussian Noise (fgn) is commonly used for that. However, if one wants to generate a long-range dependent process whose marginal law is non-Gaussian, the problem is more complex. The inverse method [START_REF] Varatkar | On-chip traffic modeling and synthesis for mpeg-2 video applications[END_REF] only guarantees an asymptotic behavior of the covariance function. We have developed, for several common laws (exponential, gamma, χ 2 , etc.), an exact method of synthesis described in [START_REF] Scherrer | Non gaussian and long memory statistical characterisations for internet traffic with anomalies[END_REF]. We can thus produce synthetic long-range dependent sample paths that can be used, as it will be shown in Section 4, in order to evaluate the impact of lrd on the performance of an on-chip network. For instance, delay sequences are not likely to have a Gaussian distribution, but rather an exponential one as we expect many small delays and few big ones. Using our synthesis method, we can produce a synthetic exponential process with long-range dependence.

Multi-phase traffic generator environment

In the following, we present our synthesis flow for building multi-phase traffic generators that can be used to replace an ip in cycle-accurate NoC performance evaluation.

Performance Evaluation Multi-Phase Traffic Generator Configuration

SoC simulation environment

We use an open source, SystemC-based, cycle-accurate and bit-accurate simulation environment: SocLib [START_REF]Soclib simulation environment[END_REF]. The environment contains cycle-accurate models for various ips. For instance, it contains a Mips R3000 processor (with its associated data and instruction cache), on-chip memories, and a component used for displaying output (referred to as tty). SocLib also includes cycle true simulators of the networks on chip developed at the LIP6 laboratory: spin and dspin [START_REF] Greiner | A generic architecture for on-chip paquets-switched interconnections[END_REF].

The application running on the Mips is composed, in addition to bootstrapping information, of the C program cross-compiled with gcc to a Mips target.

In our experiments, we use several embedded programs from the multimedia domain (still image, video and audio decoding). We make the assumption that contention only delays the communications, without changing their ordering. This is realistic for most applications and networks-on-chip.

The global simulation flow is depicted in Fig. 3. First, we generate a reference trace by running the processor ip to be emulated. This trace is obtained with an ideal network environment (no network contention), hence capturing the intrinsic communication behavior of the application. Next, we process the trace in our traffic analysis and synthesis tool explained hereafter and we obtain configuration files for our traffic generators. A generic traffic generator, further referred to as mptg, has been written once for all.

Traffic analysis

The first step of the analysis is the segmentation. It consists in identifying phases in the trace with a regular behavior. We have developed an automatic procedure for that based on existing work on cpu simulation. This part of the work is detailed in [START_REF] Scherrer | Automatic phase detection for stochastic on-chip traffic generation[END_REF]. An example of phases detected for the mp3 application is shown in Fig. 4.

On each of the identified phases, each element of a transaction: A(k), C(k), S(k), D(k) and I(k) (see 2.1), is independently analyzed with the following procedure:

• If the designer chooses deterministic traffic generation, each transaction 5uple is recorded in a file and compressed using the bz2 program, which results on average in dividing the size of the trace by 70. The resulting trace can be directly replayed by our mptg. • If the designer chooses stochastic traffic generation, a statistical analysis is performed by an automatic fitting procedure that adjusts the first and second statistical orders. The designer must choose which model they want to use (with or without lrd for instance) before the analysis can take place.

The probability distribution function (first statistical order) can be either fitted to some classical distributions (Gaussian, Exponential, Gamma, Log-Normal, etc.) or kept as it is (the model is then the probability of appearance of each value of the process). The fit is done using Maximum Likelihood Expectation and a χ 2 goodness-of-fit test is used to compare and evaluate all different solutions [START_REF] Jain | The Art of Computer Systems Performance Analysis[END_REF].

The covariance function (second statistical order) can be fitted to the one of an arma 1 process (short range correlations only), a fgn 2 (long-rangedependence only) process, or a farima 3 process (both short and long-range correlations). We use a wavelet-based estimation of the Hurst parameter [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF] widely adopted in the network traffic analysis domain (see section 2.2).

At the end of this step, a configuration file is generated. It contains, for each phase, all needed information for the traffic synthesis.

Traffic synthesis

The deterministic traffic generation is simply done by reading the compressed traffic trace.

For stochastic traffic generation purposes, we have developed an independent random number generator that can produce realizations of a wide variety of processes including non-Gaussian long-range-dependent ones (see Section 2.2.3). This generator is integrated in the mptg and the analysis simply produces the adequate mptg configuration file.

The platform designer then describes the desired platform architecture (e.g. as the one presented in Fig. 12) and uses a perl script (referred to as SocGen in Fig. 3) generating all files needed for the simulation.

At simulation time, transactions are generated by mptg according to a phase description file and a sequencer is in charge of switching between phases. At the 1 Auto-Regressive Moving Average 2 Fractional Gaussian Noise 3 Fractionnaly Intergarted ARMA end of the simulation, performance analysis indicates whether any parameters of the platform must be changed or not.

Traffic generator key features

As a summary, we highlight the main characteristics of our mptg:

• Our mptg takes into account network properties (latency, contention, etc.).

This means that the mptg will offer a good emulation of the processor independently of the network it is connected to. Hence, we can re-use the same configuration file for many different platforms.

• The global methodology presents as little manual editing as possible (e.g.

the complete simulation platform is generated). Moreover, the traffic analand synthesis flow can very easily be plugged into another synthesis environment: only the traffic generator ip needs to be adapted.

• Our mptg is able to run processes having a wide variety of parameters: deterministic replay, various stochastic models with prescribed first and second order statistics. To our knowledge this feature is not present in any cycle accurate simulators. • Our mptg is multi-phase and hence can take into account the non stationarity of the traffic.

With all these features, the designer has a very flexible tool for the design space exploration of NoC. The accuracy of the generated traffic with respect to the original trace it is suppose to emulate has been demonstrated in [START_REF] Scherrer | Generic multi-phase on-chip traffic generator[END_REF].

On any metric the mean error is less than 4%.

Experiments

In this section, we investigate the presence of long-range dependence in onchip processor traffic, and we present experimental results about the impact of long-range dependence on network-on-chip performance. In particular, we show how to use this information in order to evaluate the size of fifos in routers.

Presence of long-range dependence

In order to check for lrd in on-chip traffic, we have run several software implementations of applications from the multimedia domain: We use a simple platform (see Fig. 5) in order to truly characterize the traffic of the triplet (implementation/processor/cache). As explained in Section 3, our goal is to replace this triplet by a traffic generator. If we study communication on a more complex platform (the one in Fig. 12 for instance), the traffic of an ip is influenced by communications between other ips and by network-on-chip configuration (topology, routing protocol, etc.). Therefore, each simulation is first done on a simple platform including a Mips r3000 processor (associated with instruction and data cache), directly connected to a memory holding all necessary data represented in Fig. 5.

App. Input

Mpeg-2 2 images from a clip (176×144 color pixels)

Mp3 2 frames from a sound (44,1 kHz, 128 kbps) M-Jpeg "Lena" picture (256x256) S-Jpeg "Lena" picture (256x256) Jpeg-2000 "Lena" picture (256x256) Table 1 Inputs used in the simulations The traffic trace is recorded at the network's interface of the cache. A communication trace is extracted from the vcd (Value Change Dump) trace file, as explained in Section 3. We study the aggregated throughput, computed as the number of flits sent in consecutive time-windows of size 100 cycles. Aggregated throughput is interesting because it combines delay and size of transactions. For the sake of clarity, we present results on this quantity alone. Note that lrd is present neither in commands nor in target addresses time series.

In order to estimate the covariance function and investigate the presence of lrd in these traffic traces, we have used the wavelet-based estimator [START_REF] Abry | Wavelet analysis of long-range dependent traffic[END_REF] described in Section 2.2. Each log-scale diagram (Fig. 7-11) includes the estimated covariance (full line) and the linear regression for the estimation of the Hurst parameter (dotted line). The estimated value for this parameter, denoted Ĥ is reported in the caption.

As stated in Section 3, traffic phases can be identified, and one should pay attention, when checking for long-range dependence in a traffic trace, to consider stationary parts only. Indeed, covariance estimation tools can produce misleading results when used on non-stationary time series. Therefore, we have extracted stationary parts of the traffic traces based on information provided by our phase determination software [START_REF] Scherrer | Automatic phase detection for stochastic on-chip traffic generation[END_REF]. This is shown in Fig. 6. In each part we compute LD diagrams and estimations of the Hurst parameter as shown in Fig. 7-11. Below are the comments for each application.

• Mpeg-2 (see Fig. 7). The shape of the LD does not exhibit evidence of lrd. Indeed the estimated value for the Hurst parameter, 0.56, means no lrd is present in the trace. In this case, an iid (Independent Identically Distributed) process is a good approximation of the traffic. One can note a peak around scale 2 5 , meaning that a recurrent operation with this periodicity is present in the algorithm which might have an impact on network contention. Such behavior could be captured by an ARMA process for instance [START_REF] Brockwell | Time Series: Theory and Methods[END_REF]. • Mp3 (see Fig. 8). As for the Mpeg-2 implementation, no trace of lrd can be found in this case: the estimated Hurst parameter is close to 0.5. The other parts of the communication trace do not present any lrd either. • M-Jpeg (see Fig. 9). For this simulation, a linear behavior can be observed in the range of scales [2 4 2 14 ]. On top of this behavior, a peak is present around scale 2 11 . The behavior at larger time scales does not exhibit longrange dependence, but scale 2 14 corresponds to approximately 1.6 million cycles, which means lrd should be taken into account. For the joint modeling of lrd and small-scale behavior, farima processes can be used. • S-Jpeg (see Fig. 10). Firstly, It is interesting to note that the shape of the LD is very different from the M-Jpeg implementation, showing that the statistical properties observed are implementation rather than algorithm dependent. For the S-Jpeg implementation, long-range dependence can be observed in the limit of large scales ([2 10 2 13 ]). As for other applications, peaks are present in LD.

• Jpeg-2000 (see Fig. 11). On this application, the traffic trace exhibit a strong non-stationarity, so that the trace must be split in rather short parts for the analysis. In some of these parts, corresponding specifically to the Tier-1 entropic decoder of the Jpeg-2000 algorithm, lrd is present with an estimated Hurst parameter value between 0.85 and 0.92 (depending on the parts). In the other parts of the algorithm, no long-range dependence could be evidenced. We can conclude from these experiments that long-range dependence is not an ubiquitous property of the traffic produced by a processor associated with a cache executing a multimedia application. In some parts of the algorithms, lrd is present, however it is combined with periodicity effects which may have an equivalent impact on the network-on-chip performance.

Impact of long-range dependence

In this section we study the impact of long-range dependence on FIFO usage and communication latency in a network-on-chip. All the simulations are done on the platform shown in figure 12 in the SocLib simulation environment [START_REF] Paris | Soclib simulation environment[END_REF]. The network on chip is dspin developed at the Lip6 laboratory as an evolution of the spin network [START_REF] Greiner | A generic architecture for on-chip paquets-switched interconnections[END_REF]. It uses a mesh topology, static XY routing and wormhole memorization in the switches.

Experimental Setup

The lrd tg component is a traffic generator that produces a traffic in which the delay or inter-request time between transactions has an exponential marginal distribution and a long-range dependence of parameter H. Sizes of transactions are iid exponential random variables, with mean 10 flits (it is the typical value observed in the simulation of multimedia applications). The back tg traffic generator is used to introduce contention in the network. It injects, at a constant rate, a traffic composed of constant-sized transactions of 10 flits. The load introduced by this component is expressed as a percentage of the maximum load (1 transaction every cycle). The back tg component uses split transactions in order to truly inject a constant rate in the network. All transactions of both traffic generators are addressed to the target component: a ram (with random read or write commands) which answers in one cycle (ideal memory). The simulation length has been fixed to 10000 transactions of the lrd tg.

The major goal of this work is to be able to quantify the size of the buffers present in the switches and the wrappers (components inserted between ips and switches) as a function of the amount of lrd present in the traffic. For that, we have recorded the evolution of the usage of all fifos along the communication path used by traffic generators, as shown in figure 12. We have also computed the latency of the communications adding the latency of all fifos. We have run two types of simulations:

• fixed: in these simulations, the fifo size in the switches has been fixed to 2 flits in order to study the usage of the wrappers fifos. Indeed, with such small-sized fifos in the switches, the contention is directly transmitted to the wrappers which will be in charge of adjusting the traffic of the ip to the availability of the network. • inf: in these simulations, the fifo size in the switches and in the wrappers has been fixed to a very high (non-realistic) value of 250 flits so that the usage of all fifos (switches and wrappers) can be studied.

In both types of simulation, we have used two different communication schemes:

• no-split: in this configuration the lrd tg does not use split transaction. It waits for the reception of the k th transaction before attempting to send the (k + 1) th . The generated delay is therefore the time between the reception of the last response and the start of the new request. This communication scheme is used by peripherals and some processors with caches. • split: in this configuration the lrd tg uses split-transaction, that is to say the responses do not influence the sending of the requests. Here, the generated delay is the time between requests. This communication scheme is representative of hardware accelerators with direct memory accesses.

A particular simulation will therefore be referred to as a couple (T ,S), where T ∈ {fixed, inf} and S ∈ {no-split, split }. For instance the simulation (fixed, split) corresponds to small fifos in the switches and splittransaction communication scheme. The SoC platforms are approximately simulated at 100 000 cycles per second.

Results

Each figure contains five curves that correspond to different values of the Hurst parameter (H) of the lrd tg. The vertical axis is always a performance metric. It might be, for instance, the average latency over the whole execution or the usage (maximum or mean over the execution) of a particular fifo.

The horizontal axis corresponds to the load introduced by the back tg (in percentage of maximum load). The impact of lrd in the case of fixed is evidenced by the difference between Fig. 13 and 14. Recall that for these simulations, all the buffering occurs in the wrapper fifos because the sizes of fifos in switches are fixed to 2 flits. Fig. 13 clearly shows that lrd has no impact on performance in the case of non-split transactions because the five curves representing different values of H are identical (similar results occurs when observing the fifo of the back tg wrapper). This is typically what would appear if a processor ip (with cache) were running an application. Conversely, in the case of split transactions, the results of Fig. 14 confirm the results of [START_REF] Varatkar | On-chip traffic modeling and synthesis for mpeg-2 video applications[END_REF]: the lrd has a strong impact on performance and must be taken into account to correctly prototype NoC performance. The greater the H parameter is, the sooner network contention occurs.

The second set of figures (Fig. 15 and16) shows how we evaluate the size of the fifos needed in the switches. This must be linked with a static analysis of the path taken by the data. Fig. 12 highlights the fifos used by the data flowing from the ips to the target. We may for instance deduce that the North input fifo of router 10 will never be filled with more than 2 flits because all data coming into router 10 comes from a single source (router 11). This is confirmed in Fig. 16-(c). In these (inf,split) experiments (Fig. 15 and16) the size of the fifos in switches is fixed to a very large value (250 flits) to evaluate the maximum and usage of the fifos. The contention of the network occurs later than in the (fixed,split) experiments (comparing Fig. 15-(a) and Fig. 14-(a)) because there are more buffers on the data path. Note that the average latency on Fig. 15-(a) does not diverge except for H = 0.9. The latency diverges if the throughput given to communication is less than the injection rate of that communication. Recall that in these simulations, only the back tg component injection rate varies. The average load introduced by the lrd tg component is constant. Moreover, the round robin arbitration scheme in the router guarantees a minimum throughput to the communications of the lrd tg. Only if the fifos are small (see for instance Fig. 14-(a)), or if the injection rate of back tg is sufficiently high, will the latency diverge. In this particular simulation (inf,no-split), the average input rate of the lrd tg is not large enough to make the latency diverge. However it is interesting to note that when H is close to 1 the latency diverge as if the average rate of lrd tg was much bigger.

Looking at Fig. 15 If real time constraints are soft, the designer can decide to look at the mean usage of the fifos rather than the maximum usage. To have a more precise evaluation of the impact of setting this FIFO size, the designer can also re-run the simulation by fixing the size of router 11 input West FIFO to 50 for instance and check the impact on performance. As a summary of these experiments: lrd is more likely to have an impact on NoC performance when split-transactions communication scheme is used by Another original result presented here is a practical way to compute fifo usage and to precisely fix the size of the NoC fifos in order to reach a certain level of performance.

Related work

A recent and up to date review on works in the topic of networks-on-chip is done by Bjerregaard and Mahadevan in [START_REF] Bjerregaard | A survey of research and practices of networkon-chip[END_REF]. In this survey, traffic characterization for NoC performance evaluation is highlighted as a key issue. More specifically related to on-chip traffic generation, two main approaches have been studied: deterministic traffic generation, in which the objective is to exactly reproduce the traffic of a given ip [START_REF] Genko | A complete network-on-chip emulation framework[END_REF][START_REF] Loghi | Analyzing on-chip communication in a mpsoc environment[END_REF][START_REF] Mahadevan | A network traffic generator model for fast network-on-chip simulation[END_REF] and stochastic traffic generation which uses random sources in place of real ips [START_REF] Lahiri | Evaluation of the traffic performance characteristics of system-on-chip communication architectures[END_REF][START_REF] Pestana | Costperformance trade-offs in networks on chip: A simulation-based approach[END_REF][START_REF] Wiklund | Network on chip simulations for benchmarking[END_REF]. In the deterministic approach, Mahadevan et al. introduced, for instance, a trace compiler able to accurately reproduce the traffic of a processor [START_REF] Mahadevan | A network traffic generator model for fast network-on-chip simulation[END_REF]. As networks on chip are still in early stages of development, the main part of their performance evaluation is done using a random source. In [START_REF] Lahiri | Evaluation of the traffic performance characteristics of system-on-chip communication architectures[END_REF], stochastic processes are used for generating transaction sizes and transaction delay using several statistical laws (Poisson, Exponential, and Normal). However, none of these works proposes a fitting procedure to determine the adequate statistical models that should be used to simulate a given traffic. Thid et al. use self-similar traffic for on-chip prototyping in a recent work [START_REF] Thid | Flexible bus and noc performance analysis with configurable synthetic workloads[END_REF]. The work of Varatkar and Marculescu about long-range dependence in on-chip traffic is referenced in Section 2.2. Among the works addressing the fifo sizing, Hu et al. introduced a nonuniform NoC buffer space allocation algorithm [START_REF] Hu | ApplicationSpecific Buffer Space Allocation for NetworksonChip Router Design[END_REF]. Their algorithm was evaluated under various random traffic patterns as well as exponentially distributed traffic fitted to the communications of real applications. Chandra et al. proposed a methodology to size the fifos along a communication path [START_REF] Chandra | An interconnect channel design methodology for high perforance integrated circuits[END_REF]. They identified the various parameters that impact fifo sizing, among which the production rate, consumption rate and data burst size are the most important. Ogras et al. also pointed out in [START_REF] Ogras | Key research problems in noc design: A holistic perspective[END_REF] the fact that the NoC buffer sizing problem suffers from critical issues such as realistic traffic patterns. In this work, we have presented our proposition for having more realistic traffic patterns as well as an application of this to the buffer sizing problem.

Conclusion

We have studied both the presence and the impact of long-range dependence in the scope of networks-on-chip. Our results show firstly that long-range dependence is not an ubiquitous property of the traffic produced by on-chip processors running multimedia applications. Secondly, using cycle accurate level simulation of a complete SoC, we have shown that the low level protocol used by ips to communicate has an important impact on the possible performance degradation due to the presence of long-range dependence. In the case of nonsplit transaction, lrd is likely to have a very small impact on performance while in the case of split transaction, lrd greatly affects performance. We have also shown how the NoC designer can use our traffic generation environment to optimize the NoC for a particular SoC platform.

DFig. 1 .

 1 Fig. 1. Traffic modelling formalism

Fig. 2 .

 2 Fig. 2. Example of log-scale diagram (LD). The Hurst parameter is estimated by the slope of the dashed line (here, H = 0.83)

Fig. 3 .

 3 Fig. 3. Multi Phase Traffic Generator (mptg) Framework: Traffic analysis and synthesis flow

Fig. 4 .

 4 Fig. 4. Segmentation of the MP3 aggregated throughput into four phases

•Fig. 5 .

 5 Fig. 5. The simulation platform for lrd detection

Fig. 6 .

 6 Fig. 6. Aggregated throughput traces obtained from the execution of each implementation. The boxes mark the part of the trace on which the LD curves were obtained.

Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 78910 Fig. 7. LD of the traffic trace corresponding to the Mpeg-2 implementation. Ĥ = 0.56

Fig. 11 .

 11 Fig. 11. LD of the traffic trace corresponding to the Jpeg-2000 implementation. Ĥ = 0.89

Fig. 12 .

 12 Fig. 12. Platform used for evaluating the impact of LRD. Data paths along the switches of the NoC are represented as dash lines, FIFOs used by communications are represented inside the switches.

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Simulations FIXED and NO-SPLIT: average latency and usage of the LRD TG wrapper FIFO. (a: Latency, b: LRD TG wrapper FIFO (Mean), c: LRD TG wrapper FIFO (Max))

Fig. 15 .

 15 Fig. 15. Simulations INF and SPLIT: Maximum usage of the first three FIFOs on the path from LRD TG to TARGET (a: Latency, b: LRD TG wrapper FIFO, c: router 10 East out FIFO)

Fig. 16 .

 16 Fig. 16. Simulations INF and SPLIT: Maximum usage of the last three FIFOs on the path from LRD TG to TARGET: contention occurs in router 11 (input West FIFO). (a: router 11 West in FIFO, b: router 11 South out FIFO, c: router 10 North in FIFO)

  and 16, the designer can fix the sizes of different fifos. For instance, if the lrd tg ip has H = 0.8 and if the back tg ip emits at a rate of 40% the maximum rate, then we can see by looking at Fig.16-(a) that the input West FIFO of router 11 should have a size of 100 flits if we do not want any contention to appear (for instance, if hard real time is required).