Julien Chabé

Eric Aristidi

Aziz Ziad

Yan Fantei-Caujolle

Henri Lantéri

Christophe Giordano

Julien Borgnino

Catherine Renaud

Monitoring the atmospheric turbulence profile with high vertical resolution with the PML instrument

Keywords: Site-testing, Adaptive Optics, Interferometry

Wide-Field Adaptive Optics systems will be a key component of the future generation of Extremely Large Telescopes. All the components of such Adaptive Optics systems have to be precisely specified using the atmospheric turbulence parameters of the site, particularly the profile of the refractive index structure constant C 2 n (h). The monitor PML (Profiler of Moon Limb), for the extraction of the C 2 n (h) profile with high vertical resolution for nighttime and daytime conditions, has been developed and is now routinely exploited at the Calern Observatory (French Riviera). The PML instrument uses a differential method with two small subapertures mask through which the Moon limb or Sun edge are observed leading to a continuum of double stars allowing a scan of the whole atmosphere with high resolution in altitude. In addition, PML provides, in real time, a complete characterization of the atmospheric turbulence since it is able to measure the other turbulence parameters like the total seeing and the isoplanatic angle.

INTRODUCTION

Wide-Field Adaptive Optics (WFAO) have been proposed for the large field compensation in the perspective of the future Extremely Large Telescopes (ELTs). [START_REF] Le Louarn | New challenges for adaptive optics: extremely large telescopes[END_REF][START_REF] Hubin | Toward the adaptive optics for the 40 m class european elt[END_REF][START_REF] Ellerbroek | The tmt adaptive optics program[END_REF] The optimization of the performances of the WFAO techniques requires a precise specification of the different components of these systems. Some of these technical specifications are related to the optical parameters of the atmospheric turbulence such as seeing and isoplanatic angle. Recent studies [START_REF] Costille | Impact of the cn2 description on wide field AO performance[END_REF][START_REF] Carbillet | Anisoplanatic error evaluation and wide-field adaptive optics performance at dome c, antarctica[END_REF] show that it is necessary to have access to C 2 n (h) profiles with high vertical resolution for a better evaluation of the performance of a WFAO system.

Precise nocturnal and diurnal knowledge of the C 2 n (h) profiles (as well as wind and outer scale profiles) with high vertical resolution is also needed for Free Space Optics (FSO) for communications or for Time and Frequency Transfer in Two Way Coherent Optical Link 6 between earth and artificial satellites. Wave optics simulations to evaluate the dramatic effect of scintillation of the optical signal on fibre injection and the estimated Bit Error Rate (BER) 7, 8 must be feeded with those profiles. This helps for designing specific AO systems with higher bandwidth and capable of working over a high variability of atmospheric conditions. Real-time evaluation of accessibility of the atmospheric communication channel based on this kind of simulation will also be a need for future FSO networks. [START_REF] Horwath | Experimental verification of optical backhaul links for high-altitude platform networks: Atmospheric turbulence and downlink availability[END_REF] The optical turbulence forecasting, necessary for both astronomical observations scheduling [START_REF] Giordano | Atmospheric and seeing forecast: WRF model validation with in situ measurements at ORM[END_REF][START_REF] Giordano | Weather Research and Forecasting prevision model as a tool to search for the best sites for astronomy: application to La Palma, Canary Islands[END_REF] and optical links, is done using a numerical weather system coupled with an optical turbulence model. The calibration of these two models is then required to reach the best possible precision. The best way to do, is to use local meteorological measurements to parametrize the weather model, and local profiles (daytime and night time) of C 2 n (h) to constrain the optical turbulence model. Having a high vertical resolution C 2 n (h) profiles is important to consider, at best, phenomena within the surface layer.

Our instrument PML (Profiler of Moon Limb) has been developed for the extraction of the C 2 n (h) profile with high vertical resolution from Moon limb fluctuations. The Moon or Sun limb acts as a continuum of double stars with all possible angular separations required between two points to scan the atmosphere with a very fine resolution (see Fig. 1). The PML is the only instrument with the ability to acquire data during daylight conditions, this is an advantage for a better modeling of atmospheric turbulence. [START_REF] Giordano | Contribution of statistical site learning to improve optical turbulence forecasting[END_REF] Figure 1. Principle of PML measurement of the covariance of AA differential motion between two points on the Moon limb. Large (resp. small) angular separations give access to low (resp. high) altitude layers.

PRINCIPLE OF THE MEASUREMENT

The PML instrument measures the wavefront angle-of-arrival (AA) from the motion of the Moon limb (or Sun edge) image. PML uses the differential method of the Differential Image Motion Monitor (DIMM) [START_REF] Aristidi | A generalized differential image motion monitor[END_REF] through two subapertures of diameter D separated by a baseline B. Two images of the Moon or the Sun are produced at the focal plane of a 40cm diameter telescope. A Dove prism is used to produce a mirror image of the Moon/Sun from one of the two subapertures. Fig. 2 shows a photograph of the telescope in its dome at Calern Observatory. The pupil mask has two sub-apertures of diameter D = 6cm separated by a baseline B = 27cm. The baseline is parallel to the declination axis of the equatorial mount. Fig. 3 shows the optical scheme of the PML bench. The camera is a PCO Pixelfly CCD with a matrix of 640 × 480 pixels, working in the visible domain (the bandwidth FWHM is 320 -630 nm). The angular sampling of 0.57 arcsec/pixel is a good compromise between the spatial resolution and the field of view (365 arcsec). Every ∼ 4 minutes a 1024 images cube with an exposure time of 5 milliseconds is acquired to compute a C 2 n (h) profile.

THEORETICAL BACKGROUND

The theoretical background of PML is presented in details in. [START_REF] Chabé | Pml: a generalized monitor of atmospheric turbulence profile with high vertical resolution[END_REF] The observed Moon or Sun limb is parallel to the sub-apertures baseline and the instrument measures the differential AA through the two subapertures, the transverse covariance of the differential AA fluctuations ∆α between the two images from the two subapertures is then computed:

C ∆α (θ) = [∆α(r, θ 0 )][∆α(r, θ 0 + θ)] (1)
where ∆α(r, θ 0 ) and ∆α(r, θ+θ 0 ) represent the differential fluctuations of the lunar (or solar) limb image observed between subapertures of the PML and measured at the angular positions θ 0 and θ + θ 0 . This differential angular covariance is given by: 15

C ∆α (θ) = ¢ dh C 2 n (h) K α (B, h, θ) (2)
where the kernel

K α (B, h, θ) = 2 C α (θh)-C α (B -θh)-C α (B +θh
) is a triplet of normalized spatial covariances.

In the case of the von Kàrmàn model for a baseline , a sub-aperture diameter D (here 6cm), and a single layer at altitude h, C α ( ) is given in [START_REF] Avila | Theoretical spatiotemporal analysis of angle of arrival induced by atmospheric turbulence as observed with the grating scale monitor experiment[END_REF] by:

C α ( ) = 1.19 sec(z) ¢ df f 3 f 2 + 1 L 2 0 -11/6 × [J 0 (2π f ) + J 2 (2π f )] 2 J 1 (πDf ) πDf 2 , ( 3
)
where z is the zenithal angle, f the spatial frequency, D the sub-aperture diameter and L 0 the outer scale.

We used the assumption of a discrete turbulent profile to transform eq.2 into:

C ∆α (θ) = ∆h i C 2 n (h i ) K α (B, h i , θ), (4)
where K α (B, h i , θ) is the modified spatial covariance triplet and ∆h i is the thickness of the layer i. This equation takes a matrix form: Y = M.X where X and Y are both vectors corresponding respectively to the sampled C 2 n (h i ) and the covariance difference C ∆α (θ), the matrix M contain the modified spatial covariance triplet weighted with ∆h i : M = K α (B, h i , θ).∆h i . The values of h i ranges from 50 to 950m with a thickness ∆h i of 100m, from 1250m to 4750m with ∆h i = 500m, from 5500m to 14500m ∆h i = 1000m and from 16000m to 24000m with a thickness ∆h i of 2000m. The value for the outer scale L 0 in the kernel is taken equal to 20m. It has indeed very little influence on differential covariances. [START_REF] Ziad | First results of the pml monitor of atmospheric turbulence profile with high vertical resolution[END_REF] The C 2 n (h) profile is retrieved by solving an inverse problem on eq.4 under non-negativity constraint, via the minimization of a least squares criterion between measured and modelled differential covariances.

SITE TESTING AT THE CALERN OBSERVATORY WITH PML

The PML is part of the Calern Atmospheric Turbulence Station (CATS) located at the Calern plateau (1261m, 43.75°N, 6.88°E). The CATS station also include a Differential Image Motion Moniteur (G-DIMM) and has been designed to be autonomous and it is equipped with an allsky camera to detect the presence of clouds in the sky and a weather station for monitoring, temperature, humidity, dew point and wind velocity. All these sensors information is treated by a software devoted to observation management. While observing the PML sends its data to a website available for all the users of the observatory. The webpage shows, as a function of time, meteorological data, a picture of the night sky, the fraction of sky covered with clouds, the seeing, the scintillation and the isoplanetism angle and the C 2 n (h) profile. Over the year 2019, a total of 13126 profiles have been acquired by the PML instrument: 9224 on the Sun (during daytime) and 3902 on the Moon (nighttime).The data set is larger on the Sun, since the Moon is not always present in the sky (especially in Summer where the Moon has a low altitude above the horizon close to the full Moon). Fig. 4 shows the histogram of the number of C 2 n (h) profiles over 24 hours for the year 2019. Using the both target, Sun and Moon, the PML is able to cover almost an entire day with however a lack of data around sunrise and sunset. Fig. 5 represents the number of C 2 n (h) profiles per day for the Moon and the Sun over the year 2019. The instrument reaches a median value of 41 profiles per day for the Moon and 73 profiles for the Sun. Fig. 6 shows the probability density functions of all the layers probed by the PML for nightime (Moon) and daytime (Sun) conditions over the year 2019. We can see that PML is capable to reach an high precision and resolution in the site characterization. Several features of the C 2 n (h) profiles can be clearly identified. For example, the Sun, by heating the ground, creates an active turbulent layer is the first hundreds of meters. This induces the first layers to be localized at median value of 1.2 × 10 -14 m -2/3 during the day while its PDF decreases at 2.9 × 10 -15 m -2/3 during the night and has a larger width. The PDFs of the layers between 5km and 9.5km are also very different between day and night. During daytime, the PDFs are much larger and spread over a large span of C 2 n values (10 -16 -10 -19 m -2/3 ) while, during nighttime, the PDFs are much narrower and the range is much smaller (3 × 10 -17 -3 × 10 -18 m -2/3 ).

CONCLUSION

This paper has presented the capabilities of PML for monitoring the atmpospheric turbulence at a given site: the Calern Observatory. The PML is a C 2 n (h) profiler with high resolution in the atmosphere. This instrument is based on the measurement of the angular covariance of the AA along an extended object like the Moon or Sun. The observation statistics and some examples of the data acquired at Calern Observatory show that PML is a very good tool for turbulence monitoring and site testing. The instrument provides a detailed description of the C 2 n (h) profile evolution over the year. This will help for the scheduling of the astronomical observations over the year and gives a precious input for turbulence forecasting models.

ACKNOWLEDGMENTS

A special thanks is addressed to technical team of the Lagrange Laboratory, particularly Christophe Bailet for his precious help and to the technical staff of the Calern observatory. We thank Flavien Blary and Mohammed Marjani for their contributions at the early stages of the project. The development of the PML instrument has been supported by the French National Institute for Universe Sciences in collaboration with the Observatoire de la Côte d'Azur and the University of Nice Sophia Antipolis with financial support from CNES and Federation Doblïn. This work was supported by the LABEX Cluster of Excellence FIRST-TF (ANR-10-LABX-48-01), within the Program "Investissements d'Avenir" operated by the French National Research Agency (ANR).

Figure 2 .

 2 Figure 2. Left: the Profiler of Moon Limb instrument inside its 12ft AllSky dome. Right: The collapsible solar filters for Sun observation.

Figure 3 .

 3 Figure 3. The optical device of the PML instrument.

Figure 4 .

 4 Figure 4. Number of C 2 n (h) profiles acquired per hour over the year 2019

Figure 5 .

 5 Figure 5. Number of C 2 n (h) profiles per day over the year 2019

Figure 6 .

 6 Figure 6. Probability density function of C 2 n (h) for all the layers during day (Sun) and night (Moon) over the year 2019.

Further author information: (Send correspondence to J.C.): E-mail: julien.chabe@oca.eu, Telephone: (33) (0)4 93 40


Miguel Nicolau

Marc Schoenauer

Evolving Scale-Free Topologies using a Gene Regulatory Network Model

A novel approach to generating scale-free network topologies is introduced, based on an existing artificial Gene Regulatory Network model. From this model, different interaction networks can be extracted, based on an activation threshold. By using an Evolutionary Computation approach, the model is allowed to evolve, in order to reach specific network statistical measures. The results obtained show that, when the model uses a duplication and divergence initialisation, such as seen in nature, the resulting regulation networks not only are closer in topology to scale-free networks, but also exhibit a much higher potential for evolution.

I. INTRODUCTION

Scale-Free networks are complex networks which have a few highly connected nodes, while most nodes are poorly connected [START_REF] Barabási | Emergence of scaling in random networks[END_REF]. More precisely, in such networks, the connectivity of the nodes follows a power law: the proportion P (k) of nodes with degree k (i.e. that are connected to k other nodes) is roughly proportional to k -γ , for some positive real number γ, at least above a k given value.

Such network topology has been shown to exist in a variety of both artificial and biological systems [START_REF] Wuchty | Scale-free behavior in protein domain networks[END_REF], [START_REF] Jeong | The large-scale organization of metabolic networks[END_REF], [START_REF] Guelzim | Topological and causal structure of the yeast transcriptional regulatory network[END_REF], [START_REF] Van Noort | The yeast coexpression network has a small-world, scale-free architecture and can be explained by a simple model[END_REF], [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF], and has been widely studied because of its high resistance to random failure. Different generative models have been shown to create scale-free networks: in the original "preferential attachment" model, the network is gradually built, and new nodes attach preferentially to highly connected nodes [START_REF] Barabási | Emergence of scaling in random networks[END_REF]; however, this topology can also occur as a consequence of optimization processes, such as the wiring cost to existing software components (see [START_REF] Valverde | Scale-free networks from optimal design[END_REF] and references therein); finally, some artificial genome models, created through duplication and divergence, have been shown to generate networks with a power-law degree distribution [START_REF] Pastor-Satorras | Evolving protein interaction networks through gene duplication[END_REF], [START_REF] Kuo | Small world and scale-free network topologies in an artificial regulatory network model[END_REF]. However, all these models use rules that are not directly connected to the topology of the resulting network, and in particular do not offer an easy tuning of the statistical properties of the network they build. Using the last type of generative model -the generation of genomes through duplication and divergencethis paper investigates the possibility of designing scale-free networks with a given exponent for its power-law tail.

Genetic Regulatory Networks (GRNs) are biological interaction networks among the genes in a chromosome and the proteins they produce: each gene encodes a specific type of protein, and some of those, termed Transcription Factors, regulate (either enhance or inhibit) the expression Miguel Nicolau and Marc Schoenauer are with the INRIA Saclay -Îlede-France, TAO project-team, LRI -Bâtiment 490, Université Paris Sud, 91405 Orsay, France; email: [miguel.nicolau,marc.schoenauer]@inria.fr).

of other genes, and hence the generation of the protein those genes encode. The study of such networks of interactions provides many inter-disciplinary research opportunities, and as a result, GRNs provide an exciting and fast evolving field of research.

In order to study the characteristics of GRNs, many artificial systems have been designed, either through the modeling of biological data, or purely artificially; de Jong [START_REF] Jong | Modeling and simulation of genetic regulatory systems: A literature review[END_REF] provides a relatively recent overview of such researches.

One interesting research direction regarding the use of GRNs is the extraction and analysis (static or dynamic) of their regulation network. Previous work on the structural analysis of GRNs has provided many insights, of which the following are but a few examples. It has been shown that these networks can be grown through a process of duplication and divergence [START_REF] Wolfe | Molecular evidence for an ancient duplication of the entire yeast genome[END_REF], [START_REF] Kellis | Proof and evolutionary analysis of ancient genome duplication in the yeast saccharomyces cerevisiae[END_REF]; that they can exhibit scale-free and small-world topologies [START_REF] Bhan | A duplication growth model of gene expression networks[END_REF], [START_REF] Van Noort | The yeast coexpression network has a small-world, scale-free architecture and can be explained by a simple model[END_REF], [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF], [START_REF] Kuo | Network topology and the evolution of dynamics in an artificial regulatory network model created by whole genome duplication and divergence[END_REF]; that some specific network motifs, resembling those identified by biologists as building-blocks, are present within these artificial networks [START_REF] Wuchty | Evolutionary conservation of motif constituents in the yeast protein interaction network[END_REF], [START_REF] Milo | Superfamilies of evolved and designed networks[END_REF]; and that in response to diverse stimuli, the wiring of these networks changes over time, with a few transcription factors acting as permanent hubs, but most adapting their role as an answer to the changing environment [START_REF] Luscombe | Genomic analysis of regulatory network dynamics reveals large topological changes[END_REF].

The present work focuses on the analysis of the underlying network topologies of one artificial GRN model [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF], and of its use as a generative model for scale-free topologies. Both random genomes and genomes initialised through a duplication and divergence method are first analyzed with respect to statistical properties of the topology of the resulting interaction network. Then, the inverse problem is addressed: an Evolutionary Algorithm is used to evolve artificial GRNs so that the topology of the resulting network has some given statistical properties -more precisely, a scale-free topology with a given exponent. The results obtained show that genomes created through duplication and divergence are better suited for evolution, and generate networks exhibiting power-law tails, a clear sign of a scale-free topology.

This paper is structured as follows: Section II presents the GRN model used in the simulations, including the description and analysis of the duplication/divergence process used to initialize the genomes. Section III introduces the statistical tools used to assess the scale-free properties of the networks, along with the techniques to actually compute them. Section IV describes the experimental setup, the fitness measure and the results obtained when evolving GRNs to obtain scale-free network topologies. Finally Section V discusses those results and sketches some hints for future research directions.

II. THE GRN MODEL

A. Representation and dynamics

The artificial model described here is that proposed by W. Banzhaf [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF]. It is built over a genome, represented as a bit string, and assumes that each gene produces a single protein, with all proteins regulating all genes (including the gene that produced it).

A gene is identified within the genome by an Activator (or Promoter) site, that consists of an arbitrarily selected bit pattern: in this work, a 32 bits sequence whose last 8 bits are the pattern 01010101.

The 160 bits (5 × 32) immediately following a promoter sequence represent the gene itself, and are used to determine the protein this gene produces. This protein (like all proteins within the model) is a 32 bit sequence, resulting from a many-to-one mapping of the gene sequence: each of the 32 bits of the protein results from the application of a majority rule for each of the five sets of 32 bits taken from the 5 × 32 bits of the gene (see Fig. 1).

Upstream from the promoter site are two additional 32 bit segments, representing the enhancer and inhibitor sites: these are used for the regulation of the protein production of the associated gene.

X Y Z 0 1 0 0 0 1 1 1

bits bits

Enhancer site

Inhibitor site Promoter site Gene information 160 bits 32 bits Protein 32 bits Fig. 1. Bit string encoding of a gene. If a promoter site is found, the gene information is used to create a protein, whose quantity is regulated by the attachment of proteins to the enhancer and inhibitor sites.

The binding of proteins to the enhancer or inhibitor sites is calculated through the use of the XOR operation, which returns the degree of match as the number of bits set to one (that is, the number of complementary bits in both bit patterns). In general, a Normal distribution results from measuring the match between proteins and these sites, in a randomly generated genome [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF].

The enhancing and inhibiting signals regulating the production of protein p i are then calculated as:

e i , h i = 1 N N j=1 c j exp(β(u i,j -u i,max )) (1)
where N is the number of existing proteins, c j is the concentration of protein j, u i,j is the number of matching bits between the regulating site of gene i and protein j, u i,max is the maximum match achieved for gene i, and β is a positive scaling factor.

Given these signals, the production of protein i is calculated via the following differential equation:

dc i dt = δ(e i -h i )c i -Φ (2)
where δ is a positive scaling factor (representing a time unit), and Φ is a term that proportionally scales protein production, ensuring that i c i = 1, which results in competition between binding sites for proteins.

Note that this model simplifies some of the known characteristics of the biological regulatory process: all proteins are assumed to be Transcription Factors, that is, all proteins are used to regulate the expression of all genes: in other words, the model is a closed world. Also, the model uses only one enhancing and inhibiting site per gene. However, it captures interesting properties of actual GRNs, in particular through the genome construction technique.

B. Genome Construction

The technique of duplication and mutation proposed [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF] consists in creating a random 32 bit sequence, followed by a series of length duplications associated with a (typically low) mutation rate. It has been shown [START_REF] Wolfe | Molecular evidence for an ancient duplication of the entire yeast genome[END_REF], [START_REF] Kellis | Proof and evolutionary analysis of ancient genome duplication in the yeast saccharomyces cerevisiae[END_REF] that such evolution through genome duplication and subsequent divergence (mostly deletion) and specialisation occurs in nature.

Number of genes: An analysis of the resulting number of genes in a genome was first presented by Kuo & Banzhaf [START_REF] Kuo | Small world and scale-free network topologies in an artificial regulatory network model[END_REF]. For the sake of completeness, a similar technique has been used here to investigate the influence of the mutation rate on the number of genes per genome: 1000 genomes have been created using 14 duplication and divergence events, giving a genome length of L G = 32 × 2 14 = 524288. The resulting number of genes is shown in Fig. 2: if little or no mutation is used, a large proportion of genomes have no genes at all, but a few genomes have a large amount of genes. This was indeed to be expected: if the original random sequence contains the promoter pattern, or if it appears early in the sequence of duplications thanks to a lucky mutation, then a large number of genes will be created by the duplication process. Otherwise, little or no genes will ever exist in the genome sequence.

When the mutation rate increases, the number of genes rapidly converges towards a stable average range: with rates higher than 15%, the duplication technique becomes sufficiently randomised to roughly lead to the same number of genes per genome (around 900 here) as if using randomised genome bit-strings (or, equivalently, if using a mutation rate of 50% with the duplication/divergence process).

C. GRN Topologies

As seen before, all proteins within the model regulate the expression of all genes. The strength of this regulation is determined by the binary match between the protein pattern and the regulating sites of the destination gene (Eq. 1).

The resulting network of gene interactions can be drawn as a directed graph, with vertices connecting genes producing transcription factors to the genes they regulate [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF]. As all The x-axis plots the number of genes, the y-axis (height) the number of genomes having a particular number of genes, and the z-axis (depth) the mutation rate used. Genes were not allowed to overlap.

genes produce transcription factors, the graph of the resulting interaction network is a complete graph, where all nodes are linked together. However, because of the exponential nature of the interactions given by Eq. 1, small interactions will have almost no effect on the production of a given protein.

It is hence natural to establish a minimum matching strength (threshold) and to remove weaker regulation relationships. Moreover, by using different thresholds, different network topologies can be obtained. For instance, Fig. 3 and4 show the graphs of the same completely random genome for two slightly different thresholds (respectively 23 and 24). While almost all nodes are still connected on Fig. 3, increasing the threshold by one removes many connections, and the graph on Fig. 4 is only a small sub-graph of the previous one (nodes which become isolated are not shown, which explains the smaller number of genes). Note also how the increase of the threshold creates unconnected independent sub-graphs.

G36 G3 G41 G42 G25 G5 G49 G6 G18 G39 G48 G12 G7 G24 G26 G34 G44 G17 G8 G28 G45 G19 G9 G38 G10 G31 G53 G46 G37 G11 G50 G33 G13 G29 G32 G22 G15 G52 G40 G1 G43 G20 G21 G51 G30 G14 G2 G35
Fig. 3. Gene regulatory network for a genome of length L G = 32768, created using 10 duplication events and a mutation rate of 50%, at a threshold of 23 bits. Solid edges indicate enhancing interactions, dotted edges indicate inhibiting interactions.

G45 G8 G28 G50 G12 G26 G33 G13 G29 G32 G41 G38 G15 G18 G42 G1 G43 G21 G25 G49 G34 G30 G48
Fig. 4. Gene regulatory network for the same genome as in Fig. 3, at a threshold of 24 bits.

A completely different picture is that of genomes initialised through the duplication/divergence process described in Section II-B. Fig. 5 is an example of the topology of the interaction graph for such a genome, initialized with 1% mutation rate, using 16 as the connection threshold.

G27 G1 G2 G4 G5 G7 G10 G13 G21 G22 G23 G17 G18 G19 G26 G11 G20 G6 G8 G15 G24
Fig. 5. Gene regulatory network for a genome of length L G = 32768, created using 10 duplication events and a mutation rate of 1%, at a threshold of 16 bits.

The use of a low mutation rate results in a much shallower hierarchy of nodes, with a few master genes being connected to most of the other genes, regulating and/or being regulated by them. Varying the threshold used results in networks with similar dynamics: Fig. 6 and 7 depict the same genome, with higher connectivity thresholds (17 and 18, respectively). The presence of master genes is still clear, but their connectivity is obviously lower. Note also how some master genes disappear if the threshold parameter is increased. Another observation is that the thresholds generating "interesting" topologies for randomly created genomes are higher than those for genomes created with duplication and low mutation. This is because the latter exhibit a high degree of similarity in their bit patterns, leading to a lower value of u i,j , when applying the XOR operator (see Equation 1).

G27 G2 G4 G5 G13 G21 G23 G18 G19 G1 G11 G7 G22 G6 G8 G15 G17 G20 G24

D. Connectivity variance

In order to generalize the observations made on the graphs above, an approach similar to that of Kuo et al. [START_REF] Kuo | Network topology and the evolution of dynamics in an artificial regulatory network model created by whole genome duplication and divergence[END_REF] has been used here to analyse the relationship between the number of edges and the threshold: 100 genomes have been generated, using 14 duplication events, and the network connectivity (fraction of edges) has been computed for each threshold.

The network connectivity is defined as:

N C = #edges 2n 2 (3)
where #edges is the number of edges in the network, and n is the number of nodes, or genes (2n 2 is hence the maximum number of possible edges, as each node can be connected twice to any other node, including itself). Fig. 8 shows the connectivity as a function of the threshold, for mutation rates of 1%, 5%, 10%, and 50%. It is a clear illustration of the very different behaviors with respect to connectivity depending on the mutation rate used during the duplication/divergence process:

• A high mutation rate (or, equivalently, the completely random generation of the genome) creates a network which stays fully connected with a wide range of threshold values; then, there is a sharp transition from full connectivity to no connectivity (see also Fig. 9). Moreover, there is a very small variance between different networks.

• A low mutation rate creates a network which quickly loses full connectivity; however, its transition from full connectivity to no connectivity is much smoother than that of random networks. Moreover, there is very large variance between different networks generated with the same mutation rate. Fig. 8. Fraction of edges in a graph as compared to a fully connected network (and standard deviations), versus threshold parameter, based on samples of 100 genomes, created using 14 duplication events, and mutation rates of 1%,5%, 10%, and 50%.

III. SCALE-FREE TOPOLOGIES Even though the model used is overly simplified compared

to what is known of biological GRNs (as discussed in Section II-A), an interesting issue is to find out whether or not the resulting interaction network exhibits particular properties resembling those found in natural networks, such as being Scale-Free [START_REF] Wuchty | Scale-free behavior in protein domain networks[END_REF], [START_REF] Jeong | The large-scale organization of metabolic networks[END_REF], [START_REF] Guelzim | Topological and causal structure of the yeast transcriptional regulatory network[END_REF], [START_REF] Van Noort | The yeast coexpression network has a small-world, scale-free architecture and can be explained by a simple model[END_REF], [START_REF] Babu | Structure and evolution of transcriptional regulatory networks[END_REF], at least for certain values of the threshold used to prune the graph of connections. A characteristic feature of Scale-Free graphs is that the distribution of the degrees approximately follows some power law. But assessing such a distribution is not as obvious as it seems.

A. Measurement of Power Laws

Given a sample of some quantity, the typical method for measuring whether or not this quantity follows some power law consists in measuring whether the histogram of the sampled quantity at hand is roughly linear on logarithmic scales. A linear regression (using e.g. any Least-Squares method) can be used, and the slope of the best linear approximation will be the exponent -γ of the power-law. This method, however, has been shown to introduce systematic biases into the value of the exponent [START_REF] Goldstein | Problems with fitting to the power-law distribution[END_REF].

Another option is to work directly on the sample itself, rather than on the logarithms, and to use a non-linear curvefitting method (such as the Levenberg-Marquardt algorithm [START_REF] Levenberg | A method for the solution of certain problems in least squares[END_REF], [START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF]). In this case, however, the difficulty lies in choosing the correct parameters for the optimization method, and taking into account all points of the histogram, despite their very different orders of magnitude.

To address the limitations of the above methods, a specific method, called Maximum Likelihood Estimation (MLE), has been proposed [START_REF] Goldstein | Problems with fitting to the power-law distribution[END_REF], [START_REF] Newman | Power laws, pareto distributions and zipf's law[END_REF], [START_REF] Lee | Statistical properties of sampled networks[END_REF]. It seems to be one of the most stable methods for the approximation of the exponent of a power-law, and has been used here. The MLE method computes the exponent of the power-law as:

γ = 1 + n n k=1 ln P (k) P min -1 (4)
where P (k) is the proportion of nodes with degree k, P min the minimum value of all P (k) in the sample, and n the number of samples. Note that values of k for which P (k) = 0 are not taken into account. Furthermore, the standard error of this estimation can also be easily computed as being:

σ = √ n n k=1 ln P (k) P min -1 = γ -1 √ n (5)

B. Data Quality

Newman [START_REF] Newman | Power laws, pareto distributions and zipf's law[END_REF] observed that quite often, the tail of powerlaw distributions tends to be quite noisy, because of sampling errors: this is due to the fact that very few samples exist towards the high end of the distribution. This is certainly the case with the vertex degree distributions analysed here.

To tackle this problem, a technique known as logarithmic binning can be used [START_REF] Albert | Error and attack tolerance of complex networks[END_REF]. It smoothes the histogram by grouping the distribution data per ranges of k values with exponentially increasing sizes (e.g. 1, [START_REF] Wuchty | Scale-free behavior in protein domain networks[END_REF][START_REF] Jeong | The large-scale organization of metabolic networks[END_REF], [START_REF] Guelzim | Topological and causal structure of the yeast transcriptional regulatory network[END_REF][START_REF] Valverde | Scale-free networks from optimal design[END_REF], . . . ). This technique is also used in the present work.

C. Are GRNs Topologies Scale-Free?

Random genomes are, in terms of degree distribution, highly regular, in that their degree distribution is highly peaked; this in turn leads to potentially misleading good γ values (linear regression of 2 points is always perfect!). This can be seen in Fig. 9, which shows an example network extracted from a random genome. The vertex degree distribution is clearly Gaussian, even when plotted in a log/log graph; however, a least-squares regression gives the value γ = 1.59. Using logarithmic binning does not help: due to the proximity of all values, there are only two points left in the distribution, leading to an MLE estimation of γ = 2.219, but with a high error for the estimation (a small n leads to a high error, as seen in Equation 5). Networks extracted from initialised genomes give a completely different picture, as seen in Fig. 10. The initial distribution has a clear linear trend in a log/log scale, but is affected by noise towards the end; a least-squares regression gives the value γ = 0.767, as a result. However, by using logarithmic binning, the values towards the end are somewhat normalised, resulting in an MLE estimation of γ = 1.370.

The occurrence of misleading γ values with random genomes can be further observed in Fig. 11: the size of logarithmic bins with random genomes is much smaller, giving rise to misleading 'good' (high) γ values. Initialised genomes, on the other hand, give a wider spread of distribution sizes, with γ values typically in the range [START_REF] Barabási | Emergence of scaling in random networks[END_REF][START_REF] Wuchty | Scale-free behavior in protein domain networks[END_REF].

Though some graphs built from the artificial GRNs considered here exhibit some characteristics of scale-free networks [START_REF] Kuo | Small world and scale-free network topologies in an artificial regulatory network model[END_REF], their degree distribution is generally quite far from a true power law. Nevertheless, while random graphs, because of the poor spread of their degree distribution, seem to be difficult to modify toward more scale-free topologies, initialised ones are more promising as seed topologies for the evolution of scale-free topologies. The next section demonstrates that evolving networks created with the duplication/divergence process described is indeed possible, resulting in yet another method to construct networks with scale-free properties.

IV. EVOLUTION OF TOPOLOGIES

The objective of this section is to evolve GRN genomes, so that the resulting interaction network gets as close as possible to a target γ value, using a simple Evolutionary Algorithm.

A. The Evolutionary Algorithm

A population of bit-string genomes, such as the ones described in Section II-A is evolved using the simple bit-flip mutation as the only variation operator. The evolution is a straightforward (25+25)-ES, i.e. 25 parents give birth to 25 offspring, and the best 25 of the 50 parents+offspring become the parents of next generation. The only tricky part lies in the adaptive way to modify the mutation rate along evolution: its rate is initially set to 1% (per bit), and adapted in a way that is similar to the well-known 1/5 rule of Evolution Strategies [START_REF] Rechengerg | Evolutionsstrategie '94[END_REF]: when the rate of successful mutations is higher than 1/5 (i.e. when more than 20% mutation events result in an increase of fitness), the mutation rate is doubled; it is halved in the opposite case 1 .

In order to compare the evolvability of populations generated by the duplication/divergence method presented in Section II-B (with mutation rate 1%) and completely random populations (or, equivalently, populations built with the same method and mutation rate 50%), 50 independent runs of 50 generations have been performed with each of those initialisation procedures.

B. Fitness Function

Cohen and Havlin [START_REF] Cohen | Scale-free networks are ultrasmall[END_REF] have shown that a large proportion of networks displaying scale-free behaviour exhibit values of γ ∈ [2, 3], with some emphasis on the central value. In this work, a narrow interval around 2.5 is used, and values of γ in [2.4, 2.6] are considered ideal. The MLE method is used to compute an estimation of γ as described by Equation 4. The fitness function (to minimise) is therefore:

F (x) =    0 if 2.4 ≤ γ ≤ 2.6 2.4 -γ if γ < 2.4 γ -2.6 if γ > 2.6 + σ n (6)
The statistical error of MLE σ (see Equation 5) is added to the absolute difference to the target γ values, as an estimate of the quality of the measurement. It is divided by the number n of points in the logarithmic binned vertex degree distribution, in order to penalize even more highly "regular" distributions where only a few data points would remain after the logarithmic binning (as seen in Section III-C).

From each GRN individual, several networks are extracted, by varying the threshold value; only the threshold giving the best fitness score is kept.

C. Experimental Results

Fig. 12 shows the best fitness in the population averaged over the 50 runs, for both initialisations procedures. As can be seen from the figure, though random genomes start with a much lower fitness, they hardly improve it during evolution (there is a slight improvement of fitness, although not visible at the scale of the plot). As already mentioned, this is due to the highly regular degree distribution in random networks: this gives very few points after logarithmic binning, and leads to fakedly good γ values. But it then makes it hard to vary the number of connections between nodes, even by adjusting the threshold parameter (see Fig. 8).

On the other hand, duplication/divergence initialised genomes do start with a worse fitness (smaller γ values), but are able to evolve to much better fitness values. The widely spread degree distribution results in more data points in the vertex degree distribution after binning. Although resulting in a higher initial error when estimating γ, it also creates a larger set of potentially fit networks from each genome, by varying the threshold parameter (as per Fig. 8).

Another reason for the greater efficiency of duplication/divergence initialised genomes as initial population for evolution is their ability to improve fitness by varying the size of genomes. This is illustrated on Fig. 13, that shows that random genomes keep roughly the same size for all genomes in the population across evolution, with very small variance across runs; initialised genomes, on the other hand, vary their size much more, with a much higher variance across runs. Even though the mutation was equally likely to add or remove a gene during evolution (by creating or deleting the 01010101 promoter pattern somewhere on the genome), such operations rarely improved the fitness for random genomes, because of the small number of sample points that remained after binning for random genomes. These findings correlate well with the results already seen in Fig. 2. The difference in terms of evolution potential with regard to scale-freeness can further be seen in Fig. 14, that displays the effectiveness of mutations during typical runs: when starting from randomly created genomes, most mutation events are harmful after just a few generations; when using an initial population of genomes built through duplication/divergence, however, evolution lasts much longer, with neutral mutations starting to appear only after 24 generations, and with some beneficial mutations appearing as late as generation 46. Fig. 15 shows two examples of typical evolved networks extracted from initialised genomes, in a log/log plot, after binning. It can be seen that not all points follow a perfect line, but the distribution clearly has a power-law tail. Similar plots were obtained for most evolved networks.

However, this figure also highlights some of the drawbacks of the experimental setup used here:

• By using logarithmic binning, the resulting vertex degree distribution consists of only a few points;

• The use of logarithmic binning in the fitness function also hinders evolution somewhat, as small changes to a node created by bit-level mutation are "diluted" across the bin to which it belongs;

• The incorporation of the error measure in the fitness function makes networks with a small number of (well aligned) points in the vertex degree distribution act as local minima, from which it is very hard to escape, in spite of the penalisation of small networks in the fitness. Fig. 15. Example networks extracted from best genomes, after evolutionary process based on initialised genomes. Vertex degree distribution was logarithmically binned, and plotted on a log/log scale. The value of γ was calculated to be 2.4 on both cases: a power-law with the same exponent (straight line) is also plotted for comparison.

The first two points can be addressed by using larger genomes, thus generating larger networks (typical network sizes were around 800 -1000 nodes here). The last point is a trickier one; a potential solution could be to raise even more the influence of the number of points n in the fitness, or setting a minimum size for it.

V. CONCLUSIONS

The experiments presented in this paper demonstrate that it is possible to evolve some networks so they approach a scale-free topology with a given exponent, by optimizing a fitness measure that is directly connected to the topological property of the network, as opposed to the more classical generative methods where the scale-free property emerges from the rules that are used (or known to be used) to build the scale-free networks both in the biological and the artificial world. The long term result of such research can be to design a methodology for building artificial networks with precisely specified characteristics -motivated by known properties related to such statistical characteristics, e.g. the high resistance to random failure of scale-free networks.

The results presented in this paper also show that genomes created using the duplication and divergence method (with small mutation rate) described in the artificial GRN model proposed by W. Banzhaf [START_REF] Banzhaf | Artificial regulatory networks and genetic programming[END_REF] can be used as starting points to generate network topologies that are typical of scale-free networks. Indeed, these initialised genomes are far better suited for evolution than purely random networks, due to the larger range of degrees in the networks they encode, as well as to the wider choice of resulting networks they can provide by varying the threshold parameter that decides of the existence of an edge between nodes.

There are of course several issues that still need to be addressed with the current approach. First of all, it remains to demonstrate that different values for the exponent can be reached that way, and not only the popular value of around 2.5 that was used here. Also, the use of logarithmic binning results in a distribution with a small number of points. A possible solution to this problem can be to use overlapping bins, in order to artificially increase the size of the sample. More generally, much larger networks should be built to assess the statistical properties with more confidence. However, whereas it is not a problem to do more duplications in the initial phase of duplication/divergence, the issue when tackling larger networks will rapidly be that of CPU time: at the moment, a single evolution takes 4-6 hours of a recent Pentium computer (3.6GHz) for random networks, and 4-10 hours for duplication/divergence initialized networks, due to the higher number of threshold values that need to be checked for power-law distribution -and the main source of computational cost is the need to try several thresholds per genome. A possible solution might be to devise a heuristic in order to only evaluate promising threshold values. Another possible extension of this work would be to use localised mutations at gene encoding sections of the genomes only (or, equivalently, to remove all non-coding parts of the genome). While this will potentially increase the speed of evolution (by removing most neutral mutations), it will also remove the potential to add (or remove) genes. Though the number of genes did not vary greatly during the experiments presented here (Section IV-C), the influence of fixing the number of genes remains to be studied in more detail.

Fig. 2 .

 2 Fig.2. Histogram of the number of genes per genome, over 1000 genomes. The x-axis plots the number of genes, the y-axis (height) the number of genomes having a particular number of genes, and the z-axis (depth) the mutation rate used. Genes were not allowed to overlap.

Fig. 6 .Fig. 7 .

 67 Fig. 6. Gene regulatory network for the same genome as in Fig 5, at a threshold of 17 bits.

Fig. 9 .

 9 Fig. 9. Vertex degree distribution for the best network of a random genome, before (+) and after (vertical bars) logarithmic binning.

Fig. 10 .Fig. 11 .

 1011 Fig.10. Vertex degree distribution for the best network of an initialised genome, before (+) and after (vertical bars) logarithmic binning.

Fig. 12 .

 12 Fig.[START_REF] Kellis | Proof and evolutionary analysis of ancient genome duplication in the yeast saccharomyces cerevisiae[END_REF]. Average best fitness per generation across 50 independent runs, for random genomes and 1% duplication/divergence initialised genomes. Error bars plot the standard deviation across runs.

Fig. 13 .

 13 Fig. 13. Average population genome size across 50 independent runs, for random genomes and 1% duplication/divergence initialised genomes. Error bars plot the standard deviation.

Fig. 14 .

 14 Fig.[START_REF] Kuo | Network topology and the evolution of dynamics in an artificial regulatory network model created by whole genome duplication and divergence[END_REF]. Evolution of mutation effectiveness for random genomes (top) and 1% duplication/divergence initialised genomes (bottom), taken from a typical run. OK (resp. KO) mutations improved (resp. degraded) the fitness, while neutral mutations did not modify it.

Note that, because of the possibility of neutral mutations (especially with low mutation rates), if there were more than 50% neutral mutations, the rate was doubled in any case.

ACKNOWLEDGMENT

The authors would like to thank Wolfgang Banzhaf for his helpful suggestions. This work was supported by the Sixth European Research Framework (proposal number 034952, GENNETEC project).