Weize Quan

Kai Wang
email: kai.wang@gipsa-lab.grenoble-inp.fr

Dong-Ming Yan
email: yandongming@gmail.com

Denis Pellerin

Xiaopeng Zhang

Improving the Generalization of Colorized Image Detection with Enhanced Training of CNN

Keywords: Image forensics, natural image, colorized image, convolutional neural network, generalization, negative samples

Image colorization achieves more and more realistic results with the increasing power of recent deep learning techniques. It becomes more difficult to identify the synthetic colorized images by human eyes. In the literature, handcraftedfeature-based and convolutional neural network (CNN)-based forensic methods are proposed to distinguish between natural images (NIs) and colorized images (CIs). Although a recent CNN-based method achieves very good detection performance, an important issue (i.e., the blind detection problem) still remains and is not thoroughly studied. In this work, we focus on this challenging scenario of blind detection, i.e., no training sample is available from "unknown" colorization algorithm that we may encounter during the testing phase. This blind detection performance can be regarded as the generalization capability of a forensic detector. In this paper, we propose to first automatically construct negative samples through linear interpolation of paired natural and colorized images. Then, we progressively insert these negative samples into the original training dataset and continue to train the network. Experimental results demonstrate that our enhanced training can significantly improve the generalization performance of different CNN models.

I. INTRODUCTION

With the increasing popularity and sophistication of image editing technologies, it is now relatively easy to create edited images that are visually very plausible. For example, current advanced colorization algorithms, more or less leveraging the powerful capacity of deep neural networks, can automatically colorize a grayscale image to obtain a high-quality color image. Fig. 1 shows a group of images, the left-most one is the original color image, and the remaining three are colorized images produced by three state-of-the-art colorization algorithms (respectively with the name Ma [START_REF] Larsson | Learning representations for automatic colorization[END_REF], Mb [START_REF] Zhang | Colorful image colorization[END_REF] and Mc [START_REF] Iizuka | Let there be color!: joint end-to-end learning of global and local image priors for automatic image colorization with simultaneous classification[END_REF] from left to right), which take the grayscale version of the left-most image as input. It is indeed difficult to distinguish which images are colorized by naked human eyes. Although this technique brings convenience to people's live in fields like digital entertainment, it may also be maliciously used and potentially lead to security issues, such as confounding object recognition or scene understanding [START_REF] Guo | Fake colorized image detection[END_REF]. Therefore, distinguishing between natural images (NIs) and colorized images (CIs) has become an important research problem in image forensics.

Recently, Guo et al. [START_REF] Guo | Fake colorized image detection[END_REF] first considered and studied this new forensic problem. On the basis of the statistical differ-Fig. 1. From left to right: a natural image taken from ImageNet [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF]; three colorized images generated by the colorization method proposed in [START_REF] Larsson | Learning representations for automatic colorization[END_REF], [START_REF] Zhang | Colorful image colorization[END_REF], and [START_REF] Iizuka | Let there be color!: joint end-to-end learning of global and local image priors for automatic image colorization with simultaneous classification[END_REF], respectively. ence between NIs and CIs in the hue, saturation, dark, and bright channels, two methods, namely, histogram-based and Fisher-encoding-based, were designed to catch the forensic difference between NIs and CIs. After having obtained the feature vectors, they trained the support vector machine (SVM) classifiers to identify fake colorized images. Zhuo et al. [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF] greatly improved the detection performance using a CNNbased color image steganalyzer WISERNet (WIder SEparatethen-Reunion Network) [START_REF] Zeng | WISERNet: Wider separate-then-reunion network for steganalysis of color images[END_REF]. However, in the challenging scenario of blind detection, i.e., no training sample is available from "unknown" colorization methods that we may encounter during the testing phase of forensic detectors, the performance of both CNN-based [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF] and handcrafted-feature-based [START_REF] Guo | Fake colorized image detection[END_REF] methods in general decreases. Hereafter, we call this blind detection performance as generalization performance. Take Fig. 1 In order to cope with the challenging scenario of blind detection, in this paper we introduce a simple yet effective method. We construct negative samples via linear interpolation of paired natural and colorized images available in the training dataset, and iteratively add them into the original training dataset for additional and enhanced CNN training. This procedure is fully automatic, and can allow us to obtain stable and high generalization performance of the CNN.

The rest of this paper is organized as follows. Section II presents the technical details of the proposed method. Sec-tion III reports the performance evaluations for our method. Section IV draws the conclusions and proposes some future working directions.

II. PROPOSED FRAMEWORK

A. Motivation

To the best of our knowledge, there is no existing work that considers the generalization capability yet for CNN-based image forensics. In fact, this is a highly challenging scenario because no training samples of the "unknown" colorization algorithms are available. In other words, we want the trained network to be able to successfully detect colorized images generated by new colorization methods that remain unknown during the training of CNN. This is a very realistic situation which can be commonly encountered after deploying a forensic detector in practical applications. We solve this challenging generalization problem through a simple yet effective approach, i.e., inserting additional negative samples that are automatically constructed from available training samples, in order to carry out an enhanced training of CNN and thus to obtain an appropriate decision boundary for this classification problem. Besides considering the CNN model proposed in the very recent work of [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF], in this paper, we also construct a different CNN model so as to validate and show that our enhanced training can work well on different networks.

B. Architecture of Networks

In this subsection, we describe the architecture of considered networks. The first layer (with so-called SRM, Spatial Rich Model, kernels [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF]) of WISERNet is untrainable, while all weights of our designed network are trainable and thus we call it as AutoNet (Automatic Network). Let Ck(M or A) denote a Convolution-BatchNorm-ReLU(-MaxPool or -AveragePool) layer with k filters. Fk(R) denotes a fully-connected layer with k neurons (and with ReLU). The architecture of Au-toNet is C32-C64M-C128M-C256M-C256M-C512M-C512M-C512-F2. All convolutional kernel sizes in AutoNet are 3 × 3.

For layers 1-7, each convolutional layer (conv) is with the zero-padding of 1, and all max-pooling layers in AutoNet have the same kernel size of 3 × 3 and a stride of 2. For conv1, we use TanH as activation. Here we also give the architecture of WISERNet and more details can be found in [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF]. The architecture of WISERNet is SRM-C72A-C288A-C1152A-F800R-F400R-F200R-F2, where SRM refers to channel-wise convolution where the convolutional kernels are fixed as the thirty 5 × 5 SRM filters borrowed from [START_REF] Fridrich | Rich models for steganalysis of digital images[END_REF].

C. Negative Sample Insertion

According to our observation, there is a certain degree of performance decrease in the challenging blind detection scenario, not only for traditional handcrafted-feature-based methods [START_REF] Guo | Fake colorized image detection[END_REF], but also for CNN-based approaches (AutoNet and WISERNet [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF]), although the latter has better performance. In details, for a traditional or CNN-based model trained on dataset constructed by one specific colorization algorithm, the test performance on datasets constructed by other colorization algorithms is sometimes limited for colorized images. The possible reason of this performance drop is that colorized images produced by a specific colorization algorithm tend to be equipped with a particular internal property, but CIs of different colorization algorithms are very likely to have different properties.

To clearly illustrate the encountered problem with an example, we train the AutoNet on the dataset constructed by colorization method Mb [START_REF] Zhang | Colorful image colorization[END_REF], and test on the datasets constructed by Ma [START_REF] Larsson | Learning representations for automatic colorization[END_REF] and Mc [START_REF] Iizuka | Let there be color!: joint end-to-end learning of global and local image priors for automatic image colorization with simultaneous classification[END_REF], respectively. It should be noted that Ma and Mc are the "unknown" colorization algorithms, and thus the corresponding samples of Ma and Mc are not used in the training process. We use t-distributed stochastic neighbor embedding (t-SNE) [START_REF] Van Der Maaten | Visualizing high-dimensional data using t-SNE[END_REF] to project the highdimensional deep features (the output of conv8 of AutoNet, and its dimension is 512) of testing data constructed by above three colorization methods onto the two-dimensional map, and detailed visualization results are shown in Fig. 2. Comparing Fig. 2(a), (b) and (c), we find that the distributions of NIs (red squares) are relatively stable with a rather high intra-class variation, which is somehow expected; in the meanwhile, CIs (blue symbols) are more tightly clustered for each colorization algorithm but their locations change a lot for different methods [please compare the CIs in (a), (b) and (c), which correspond to Mb, Ma and Mc, respectively]. This is reasonable because the different colorization methods tend to have not exactly the same internal characteristics and hence the corresponding CIs have different locations in the feature space. When the features of CIs produced by "unknown" colorization algorithms (here Ma and Mc whose samples are not used for training) are near the decision boundary of the CNN (which is trained by using NIs and CIs produced by a "known" colorization algorithm, here Mb), and at the same time the decision boundary is relatively close to colorized images, there are high probabilities to misclassify the "unknown" CIs. For instance, many CIs in Fig. 2(b) (blue circles with red + in the figure) are wrongly predicted as NIs.

We would like to find a simple yet effective method to solve the encountered problem. The idea is that we make use of the available training samples (and only these samples) to construct an appropriate decision boundary which can lead to better generalization performance. A feasible and intuitive solution is to add negative samples (with same labels as CIs) near the initial decision boundary of the CNN, so as to make the CNN be more "strict" about the predictions of CIs and somehow push the classification boundary towards NIs. As such, it is expected that the "unknown" CIs located close to the initial decision boundary [e.g., those shown in Fig. 2(b)] have more chance to be correctly classified with the new classification boundary which would be closer to NIs. More precisely, we construct negative sample through linear interpolation between paired NI and CI which share the same grayscale version and only differ in chrominance components. The corresponding formulation is shown below: As analyzed above, adding negative samples and conducting additional training will push the classification boundary towards NIs. Thus, the classification accuracy on the NIs will gradually decrease as more and more negative samples are inserted. The classification accuracy of network on validation dataset also slightly decreases because the "known" CIs are almost all correctly classified and this accuracy mainly depends on the classification accuracy on the NIs. However, in the meanwhile the CIs constructed by "unknown" colorization algorithms are expected to be classified more correctly, implying a better generalization capability. Obviously, there is a trade-off between the classification accuracy (on data similar to the training samples) and generalization performance (mainly on "unknown" CIs) for the network. Therefore, without being able to directly measure the generalization during training of network, we consider the classification accuracy on NIs (on the so-called natural validation dataset V) as a measure to select the final model in the process of additional training with negative sample insertion. In our work, we design a thresholdbased model selection criterion. This threshold (θ) essentially determines the degree of final classification accuracy that can be accepted by user or current task. Generally speaking, larger θ means that the selected model has less high classification accuracy, but better generalization performance. Basically, we set θ = β • error rate, where β is a user defined parameter and error rate is the classification error rate (in %, measured on natural validation dataset V) of the CNN model trained with the original training dataset D before negative sample insertion. This criterion simply defines the maximum tolerable value of the relative increase of error rate on V induced by enhanced training. In our experiments, we set β = 2. One exception is that when error rate is very small (less than 1%), we set θ = 2%, meaning that we can slightly relax the constraint on classification error rate to obtain relatively large improvement of generalization performance.

I N S = α • I N + (1 -α) • I C , (1)
Algorithm 1 illustrates the training process with negative sample insertion. It is worth noting that we only use CIs of a "known" colorization method but in a better way to construct a more appropriate decision boundary. In our experiments, this insertion is an iterative process with four iterations, i.e., the α is increased from 0.1 to 0.4 with step of 0.1. Given a CNN model M trained by using original dataset D, and some basic construct negative samples from P using Eq. ( 1) and insert them into N . if INS is misclassified then

9:
remove corresponding pair from P. settings for CNN training, such as initial learning rate lr 0 and S epochs for each insertion, we first compute error rate on V and then the threshold θ, which is used for final model selection. For each round of negative sample insertion, we construct negative samples and insert them into the dataset D. Then, we update the parameters of model M using new training dataset, and compute the error rate on V starting from the second half of training process (i.e., from S 2 -th epoch for each insertion, where . is the integer ceiling operator), because from that time the model becomes relatively stable. After each insertion, we test the negative samples produced by previous iteration. If a negative sample is misclassified, i.e., the predicted label is NI and not consistent with its ground-truth label, then we stop using the corresponding pair to construct negative sample (i.e., we remove corresponding pair from P as described in line 9 of Algorithm 1). In fact, this operation can slightly reduce the amount of negative samples, and does not weaken the performance of the network. After four iterations of insertion, we select the final CNN model. It is worth mentioning that when α 0.5, the negative samples will be close to NIs, and this is likely to have more impact on the classification of NIs. Here we take a conservative and experimentally effective approach, i.e., stopping the negative sample insertion process after four iterations.

The complete training process of CNN model includes two stages: (1) using the original training dataset to train the deep model from scratch until convergence (normal training); (2) iteratively adding new negative samples into the original training dataset and continuing to train the model as summarized in Algorithm 1 (enhanced training). Fig. 3 shows the error rate curves of a complete training process of AutoNet. In the first

III. EXPERIMENTAL RESULTS

A. Implementation Details

All the experiments are implemented with PyTorch 0.3.1 [10]. The GPU version is GeForce ® GTX 1080Ti of NVIDIA ® corporation. All images in our experiments are resized to 256 × 256 using bicubic interpolation, and for each image, we convert its pixel values to [-1, 1] (we first rescale the pixel values from the range [0, 255] to the range [0.0, 1.0], and then subtract these values by 0.5 and divide by 0.5). Stochastic gradient descent (SGD) with a minibatch of 20 is used to train AutoNet. Each minibatch contains 10 natural images and 10 colorized images. We randomly shuffle the order of training dataset after each epoch. For SGD optimizer, the momentum is 0.9 and the weight decay is 1e-4. The base learning rate is initialized to 1e-4. For the normal training (only using original training dataset) of AutoNet, we divide the learning rate by 10 every 20 epochs, and the training procedure stops after 60 epochs. For the normal training of WISERNet, we follow the setting described in [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF]. As shown in line 13 of Algorithm 1, for the enhanced training of AutoNet and WISERNet, we adopt the same strategy about learning rate: the learning rate is divided by 10 every 15 epochs (it is enough to guarantee the convergence after new negative sample insertion), and the training procedure stops after 60 epochs, i.e., 4 iterations of negative sample insertion.

Following [START_REF] Guo | Fake colorized image detection[END_REF] and [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF], we also employ the half total error rate (HTER) to evaluate the performance of the proposed method. The HTER is defined as the average of misclassification rates (in %) of NIs and CIs. In this work, all reported results are the average of 7 runs.

B. Effect of Negative Sample Insertion

Before evaluating the proposed method, we provide the details of datasets used in our experiments. Following [START_REF] Guo | Fake colorized image detection[END_REF] and [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF], three state-of-the-art colorization algorithms, Ma [START_REF] Larsson | Learning representations for automatic colorization[END_REF], Mb [START_REF] Zhang | Colorful image colorization[END_REF], and Mc [START_REF] Iizuka | Let there be color!: joint end-to-end learning of global and local image priors for automatic image colorization with simultaneous classification[END_REF] are adopted for producing CIs. NIs come from ImageNet dataset [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF]. We use 10,000 natural images from ImageNet validation dataset to construct training dataset and validation dataset, and the ratio is 4:1. The exact indexes of these images are shared by the authors of [START_REF] Larsson | Learning representations for automatic colorization[END_REF]. Then, we remove the 899 grayscale images and 1 CMYK (cyan, magenta, yellow, and black) image from the remaining 40,000 images of ImageNet validation dataset (the total number of images in this dataset is 50,000), and obtain 39,100 natural images to construct testing dataset. Note that, the magnitude of testing dataset is far larger than the settings reported in [START_REF] Guo | Fake colorized image detection[END_REF] and [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF]. We employ the three colorization methods mentioned above to produce the corresponding colorized images.

In this paper, we propose negative sample insertion to improve the generalization performance of CNN-based detectors. As described in Section II-C, this enhanced training uses natural validation dataset V to select the final model, and we randomly select 20,000 NIs from ImageNet test dataset [START_REF] Deng | ImageNet: A large-scale hierarchical image database[END_REF] to construct V. Table I reports the performance of AutoNet and WISERNet before (i.e., the rows of "AutoNet" and "WISER-Net") and after (i.e., the rows of "AutoNet-i" and "WISERNeti") negative sample insertion. We do not present the results of handcrafted-feature-based methods proposed in [START_REF] Guo | Fake colorized image detection[END_REF] because as shown in [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF] and also verified by our experiments, CNNbased method has significantly better performance in terms of both accuracy and generalization. The difference between the results of the row of "WISERNet" and those reported in [START_REF] Zhuo | Fake colorized image detection with channel-wise convolution based deep-learning framework[END_REF] is probably due to the differences in the generation of experimental data and the number of testing images (we use much more testing data). It is worth mentioning that here we focus on the generalization improvement after applying our proposed enhanced training for the two networks (i.e., Au-toNet and WISERNet), rather than the performance difference between them. We leave the architecture comparison and the design of CNN of better generalization as a future work. From Table I, we can see that the effect of negative sample insertion, i.e., improving the generalization of network, is consistently stable for these two networks (except for one case, trained on Mc and tested on Mb for WISERNet, but with a very low final error rate of 1.08%). The negative sample insertion leads to slight decrease of the classification accuracy, however, the generalization performance of network usually has apparent improvement. For example, the initial generalization error of WISERNet trained on Mb and tested on Ma is 33.30%, and then reduces to 4.74% after enhanced training using negative samples, with a slight increase of classification error from 0.16% to 0.94%. This is also consistent with previous analysis (Section II-C) that there is a compromise between the accuracy and the generalization performance, and our negative sample insertion method can achieve a satisfying trade-off. Very lately, we became aware of a recently proposed "mixup" learning principle [START_REF] Zhang | mixup: Beyond empirical risk minimization[END_REF] which regularizes the neural network and encourages the trained model to behave linearly in-between training examples. Although the linear interpolation is also used, there is an essential difference: "mixup" results in the linearly-transitioned decision boundary, while our method pushes the decision boundary towards NI. Based on the respective standing point, for the linear interpolation itself, [START_REF] Zhang | mixup: Beyond empirical risk minimization[END_REF] uses the interpolation factor in the range of [0, 1] to combine pair of raw inputs and their labels, whereas our method uses that of {0.1, 0.2, 0.3, 0.4} (forensically negative) and sets the label of new generated image as CI (the socalled negative sample). In addition, "mixup" is a form of data augmentation that implicitly affects the generalization of network, whereas our enhanced training explicitly controls the decision boundary and then improves the generalization of CNN-based detectors. In order to compare the "mixup" and our method, we train the model with "mixup" where the learning rate schedule is exactly the same as the normal training of AutoNet and the results are shown in Table I (the row of "AutoNet-mixup"). We set the "mixup" hyperparameter α = 0.4 as recommended in [START_REF] Zhang | mixup: Beyond empirical risk minimization[END_REF]. Obviously, the generalization of our enhanced training based on negative sample insertion is significantly better than that of "mixup", only with a slight decrease of the classification performance (please compare the rows of "AutoNet-i" and "AutoNet-mixup").

At last, we visualize deep features of AutoNet-i using t-SNE [START_REF] Van Der Maaten | Visualizing high-dimensional data using t-SNE[END_REF], and the results are shown in Fig. 4. Here, deep features are the output of conv8 of AutoNet-i, and its dimension is 512. The corresponding visualizations of the model before negative sample insertion are shown in Fig. 2. The testing data is also the same in Fig. 4 and Fig. 2. By comparing the border of correctly classified CIs, i.e., blue symbols with a blue + inside, in Fig. 2(d) and Fig. 4(d), we can find that the latter has fewer misclassified CIs, and This confirms that our negative sample insertion scheme can push the decision boundary towards NIs to some extent and accordingly improve the generalization performance.

IV. CONCLUDING REMARKS

In this paper, we considered the challenging blind detection scenario and proposed an effective method based on negative sample insertion to improve the generalization capability of CNN-based models. The generalization performance is noticeably and consistently improved, with a very slight decrease of the classification accuracy. Our source code is available at https://github.com/weizequan/NIvsCI.

In the future, we are interested in employing the proposed enhanced training to improve the generalization performance of other kinds of forensic methods whenever applicable. We would like to explore other approaches to understanding and enhancing the generalization capability of neural networks.

  as an example, the second and fourth images (produced by Ma and Mc) are misclassified as NI by a CNN model trained on NIs and CIs generated by Mb. In the meanwhile, although not being very rigorous, we choose to use the term "classification accuracy/performance" to indicate the detection performance on testing data in which CIs are generated by a same colorization method known by the training procedure.

Fig. 2 .

 2 Fig. 2. The deep feature visualization with t-SNE [9]. The model is trained on the original dataset where CIs are generated by Mb. "C" means colorized images and "N" means natural images. "C-X" means the colorized images produced by X colorization method, for example, "C-Ma" corresponds to CIs generated by Ma colorization algorithm. "Y-pred" means that the predicted label of CNN is Y. We randomly select 900 natural images from validation dataset splitting them into three equal subsets of 300 images, and then we construct corresponding colorized images using Mb, Ma, and Mc for every 300 images. The deep feature is the output of conv8 of AutoNet, and the dimension is 512. (d) is the combination of (a) [Mb], (b) [Ma], and (c) [Mc].

Algorithm 1

 1 Enhanced training of CNN model with negative sample insertion Input: M, lr 0 , S, V, D and the set of corresponding natural and colorized image pairs P constructed from D. Output: final model after enhanced training. Initialization: current learning rate lr = lr 0 , set of negative samples N = ∅, set of error rates on V of candidate CNN models R = ∅. 1: compute error rate of M. 2: compute θ. 3: for all α ∈ {0.1, 0.2, 0.3, 0.4} do 4:

5 :

 5 update training dataset: D = D ∪ N . 6:update the parameters of M for S epochs. In the second half of training process, compute error rate on V for each model, and insert this value at the end of R.

7 :

 7 for all INS ∈ N do 8:

  current learning rate: lr = lr • 0.1. 14: end for 15: select i-th model which satisfies max i {ri|ri ∈ R, ri < θ}.

Fig. 3 .

 3 Fig. 3. Error rate curves of a complete training of AutoNet. The network is trained on Mb [2], and tested on Ma [1] and Mc [3]. The error rates (in %) on CIs produced by these three methods are shown in (b), (c), and (d), respectively. The error rate on V is shown in (a). Black dotted line separates the two stages of normal training (60 epochs) and enhanced training (4×15=60 epochs). The green circle in (a) stands for the final selected model.

Fig. 4 .

 4 Fig.4. The deep feature visualization of AutoNet-i with t-SNE[START_REF] Van Der Maaten | Visualizing high-dimensional data using t-SNE[END_REF]. The model is obtained through enhanced training of the previously trained model (used in Fig.2). The meaning of symbols is same as that of Fig.2. It is worth noting that in t-SNE the transformation used for dimension reduction and the obtained visualization depend on the input data. Therefore, transformation and visualization in this figure are different from those of Fig.2.

TABLE I THE

 I PERFORMANCE (HTER, IN %, LOWER IS BETTER) OF THE TWO CNN-BASED METHODS (AUTONET AND WISERNET [6]) ON IMAGENET VALIDATION DATASET [5]. FOR THE SAKE OF CLARITY, THE GENERALIZATION PERFORMANCE RESULTS ARE PRESENTED IN ITALICS.

	Method	Ma	Mb	Mc
		Ma Mb Mc Ma Mb Mc Ma Mb Mc
	AutoNet	0.56 10.57 10.62 31.65 0.19 6.16 13.93 1.91 0.72
	AutoNet-i 1.02 6.94 5.12 5.13 0.94 1.92 3.33 1.75 1.14
	AutoNet-mixup 0.89 12.45 15.35 20.68 0.34 10.04 8.42 2.25 0.76
	WISERNet 0.29 2.21 10.74 33.30 0.16 7.88 5.80 0.59 0.36
	WISERNet-i 0.98 1.22 2.29 4.74 0.94 2.04 2.46 1.08 0.98

ACKNOWLEDGMENT

This work is funded in part by Beijing Natural Science Foundation (L182059), National Natural Science Foundation of China (61620106003 and 61772523), and the French National Research Agency (DEFALS ANR-16-DEFA-0003, ANR-15-IDEX-02). W. Quan acknowledges the support from the UCAS Joint PhD Training Program. We would like to thank Dr. G. Larsson, Dr. R. Zhang and Dr. S. Iizuka for kindly sharing the source code of their colorization algorithms described respectively in [1], [2] and [3], Dr. L. van der Maaten for making the t-SNE tool [9] publicly available, L. Zhuo for the training details of their network [6] and Dr. Y. Guo for detailed and helpful discussions about their work [4].


B Girardot

M Gauduchon

F Ménard

J C Poggiale

Does evolution design robust food webs ?

Keywords: community evolution models, environmental richness, robustness, connectance

Theoretical works that use a dynamical approach to study the ability of ecological communities to resist perturbations are largely based on randomly generated ecosystem structures. In contrast, we propose here to asses the robustness of food webs drawn from ecological and evolutionary processes with the use of community evolution models. In a first part, with the use of Adaptive Dynamics theoretical framework, we generate a variety of diversified food webs by solely sampling different richness levels of the environment as a control parameter, and obtain networks that satisfactory compare with empirical data. This allows us to highlight the complex, structuring role of the environmental richness during the evolutionary emergence of food webs. In a second part, we study the short-term ecological responses of food webs to swift changes in their customary environmental richness condition. We reveal a strong link between the environmental conditions that attended food webs evolutionary constructions and their robustness to environmental perturbations. When focusing on emergent properties of our evolved food webs, especially connectance, we highlight results that seem to contradict the current paradigm. Among these food webs, the most connected appear to be the less robust to sudden depletion of the environmental richness that constituted their evolutionary environment. Otherwise, we appraise the "adaptation" of food webs, by examining how they perform after being suddently immersed in an environment of modified richness level, in comparison with a trophic network that experienced this latter environmental condition all along its evolution.

Introduction

Understanding and linking both the structure and functioning of ecosystems has been addressed as one of the major issues in ecology, particularly due to increasing pressures on biodiversity [START_REF] May | Unanswered questions in ecology[END_REF][START_REF] Thompson | Food webs: Reconciling the structure and function of biodiversity[END_REF]. The large diversity/complexity-stability debate (see [START_REF] Shear | The diversity-stability debate[END_REF][START_REF] Landi | Complexity and stability of ecological networks: A review of the theory[END_REF] for reviews), rooted in this issue, has been approached from different angles. A first angle tackled this complexity by studying food webs -the representation of who eat whom in ecosystems -since they provide tractable depictions of biodiversity and ecosystem structure and function [START_REF] Dunne | Network structure and biodiversity loss in food webs: Robustness increases with connectance[END_REF]. Pioneer food web models such as the cascade model [START_REF] Joel E Cohen | Community food webs: data and theory[END_REF] or the niche model [START_REF] Williams | Simple rules yield complex food webs[END_REF] were able to produce satisfactory food web structures, with the ability to transpose some network properties from natural food webs. However, among the limitations of these models is first the fact that they do not propose an explanation of how these properties emerge, and second that they do not provide a quantification of interactions strengths, prerequisite for a population dynamics model of interacting species. Past decade has seen the development of community evolution models, that came to adress these limitations [START_REF] Loeuille | Evolutionary emergence of size-structured food webs[END_REF][START_REF] Ingram | Niche Evolution, Trophic Structure, and Species Turnover in Model Food Webs[END_REF][START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF]. They provide demographic and evolutionary dynamics, presenting resulting food webs as the "winnowed survivors of evolutionary processes" [START_REF] May | Unanswered questions in ecology[END_REF]. In those models, species are characterized by some traits (unique to multiple, depending on studies), which drive the demography of the whole interacting community on a given time scale. On another time scale, these traits evolve according to the selection pressures exerted on each of them. New species are introduced randomly, ultimately leading to the emergence of communities. The diversity and structure of the resulting food webs are hence emergent properties controlled by small scale species interaction parameters. The model developed by Loeuille and Loreau (8) has the advantage to be relatively simple as it implies only one characterizing trait for a species, the body size at maturity. However, networks that emerge from this model have been shown to underlie a very uniform regular structure, in particular because species have the same feeding characteristics (same feeding distance and same specialization level). In order to reproduce a more realistic diversity in feeding strategies and hence in food webs structure, it is necessary to consider these feeding parameters as evolving traits. Several studies do have attempted to consider the further evolution of these traits (e.g. 11). Still with the aim of explaining the structure and functioning of food webs, it was noted that processes that operate at global scale (e.g environmental productivity, climate) have rarely been shown to be determinant for the structure of communities, even if some exceptions exist [START_REF] José | 17 integrating species interaction networks and biogeography[END_REF]. Yet, several hypotheses have been formulated to explain how community structure could be related to environmental attributes, such as productivity and/or ecosystem sizes [START_REF] Post | Ecosystem size determines foodchain length in lakes[END_REF]. Among properties of community structure, food chain length, or connectance, has been found to vary along environmental gradients, in lakes [START_REF] Post | Ecosystem size determines foodchain length in lakes[END_REF], streams [START_REF] Thompson | Energy availability, spatial heterogeneity and ecosystem size predict food-web structure in streams[END_REF] or high latitude marine ecosystem in a more recent study [START_REF] Kortsch | Food-web structure varies along environmental gradients in a highlatitude marine ecosystem[END_REF]. To Girardot | bioRχiv | June 6, 2019 | 1-10 our knowledge, no studies using community evolution models have focused on the effect of parameters influencing these global-scale processes on the emergent structures, either because it was not their main topic, or because they led only to expected results. This is the case, for example in the seminal study of Loeuille and Loreau (8), in which discussions on the effects of variations in a parameter affecting the productivity of the environment in which the emerging network evolves (i.e global scale parameter) are brief, because they result in monotonous trends: an increase of specific diversity and maximum trophic level with richness of the evolutionary environment. Among the subsequent, more detailed community evolution model, such as Allhof's model [START_REF] Korinna | When do evolutionary food web models generate complex networks[END_REF], none has investigated so far the effects of these global scale parameters on the structure of resulting communities. Although they are still in their infancy, community evolution models show great potential, especially because they can bridge the gap between ecosystem and community ecology [START_REF] Loreau | Linking biodiversity and ecosystems: Towards a unifying ecological theory[END_REF]. Among the promising but still unexplored avenues in the exploitation of such models is the analysis of both short-term ecological responses and long-term evolutionary responses of communities to disturbances [START_REF] Loreau | Linking biodiversity and ecosystems: Towards a unifying ecological theory[END_REF].

A last angle mentioned here to broach the diversitystability debate is to focus on how to relate the structure of food webs to their fragility. For example, motivated by the observations of increasing species losses [START_REF] Barnosky | Has the Earth's sixth mass extinction already arrived?[END_REF][START_REF]Millennium Ecosystem Assessment. Millennium ecosystem assessment (ma)[END_REF], several studies have investigated the links between primary species losses on communities, the so called secondary extinctions, and the structures of food webs [START_REF] Dunne | Network structure and biodiversity loss in food webs: Robustness increases with connectance[END_REF][START_REF] Ebenman | COMMUNITY VIABILITY ANALYSIS: THE RESPONSE OF ECOLOGICAL COMMUNITIES TO SPECIES LOSS[END_REF][START_REF] Solé | Complexity and fragility in ecological networks[END_REF][START_REF] Eklöf | Species loss and secondary extinctions in simple and complex model communities[END_REF].

Here, we try to establish the link between these different approaches to the diversity-stability debate. In the first part, we build up a model inspired from [START_REF] Loeuille | Evolutionary emergence of size-structured food webs[END_REF], but modified by including two additional adaptive traits in order to bring out diversified feeding behaviors. The evolution component of the model is also improved and take opportunity of the Adaptive Dynamics toolbox. We then interested ourselves in the effects of the variation of a single control parameter that determines the richness of the environment during evolution (global scale processes parameter) and its effects on emerging structures. The properties of the food webs obtained in this way are then confronted to a large set of empirical data [START_REF] Cohen | Ecoweb 1.1: Ecologists' cooperative web bank[END_REF]. In the second and main part of the article, we study the responses on the ecological timescale of the evolved food webs to abrupt changes in the richness of their evolutionary environment. We quantify the effect of the level of environment richness that constituted the evolutionary environment of a food on its robustness when it is subsequently confronted to a different condition. We then assess and discuss the possible direct links between emerged structure properties of evolved networks and their robustness. Lastly, we investigate whether an evolved food web performs the best in its evolution-ary condition, compared to other networks that have first evolved in different richness levels and abruptly undergone this former condition. It contributes to the debate of considering the analogy of a food webs as a super-organism, that would become the fittest to its environment.

Model

Trophic network. The model is derived from Loeuille and Loreau's work [START_REF] Loeuille | Evolutionary emergence of size-structured food webs[END_REF] (LL in the following), with, as we shall see later on, some modifications both in the ecological hypotheses and in the methods we used. For species i in a community, we denote x i its typical (average, maturity) body size. Its feeding kernel is Gaussian-shaped and depends on two parameters: the preferred body size difference d i between its own body size and the size of the species it consumes, and s i the width of the kernel, so that it optimally consumes species of body size x id i , with s i determining its degree of generalism as a consumer. This setting allows for the determination of the connected species in the trophic network but also provides a weight on the connections. Thus, interactions in the trophic network are not set a priori but emerge from the relative body sizes of every species pair in the community, combined with their feeding kernels. As in subsequent papers (e.g. 10), we removed some limitations present in the original model.

Adaptive Dynamics.

For each species, we let the body size x i and also the feeding kernel parameters d i and s i be adaptive traits that undergo evolution. We here depart from LL in which only the body size evolves, following some other authors that modified the original model [START_REF] Ingram | Niche Evolution, Trophic Structure, and Species Turnover in Model Food Webs[END_REF][START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF][START_REF] Korinna | When do evolutionary food web models generate complex networks[END_REF][START_REF] Kt Allhoff | Biodiversity and ecosystem functioning in evolving food webs[END_REF]. The evolution of traits is driven by the selection pressures, which depend in turn on the network structure. But since network structure is determined by species traits, it changes as well on the evolutionary time scale and in turn reshapes the selection pressures on species. This is a school case of coevolution [START_REF] John N Thompson | The coevolutionary process[END_REF]. Adaptive Dynamics theoretical framework [START_REF] Metz | How should we define 'fitness' for general ecological scenarios?[END_REF][START_REF] Dieckmann | The dynamical theory of coevolution: a derivation from stochastic ecological processes[END_REF][START_REF] Stefan | Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree[END_REF] is particularly well suited to encompass such situations in which co-evolution of traits can only be understood from a precise description of the selection pressures and where it feeds back on them. Adaptive Dynamics provides us a very powerful tool that allows us to model efficiently the coevolution of adaptive traits through a set of ordinary differential equations for expectation of traits values: the canonical equation.

Population dynamics.

Adaptive Dynamics key ingredient is the supplying of the invasion fitness fonction that allows to assess the invasion success of a mutant, according to its new trait value and the state of the resident community in which it appears. Invasion fitness function is built up on the explicit population dynamics of interacting species that both allows to determine the steady state of populations in the community and the initial dynamics of the mutant lineage just after it appears. LL provides such population dynamics model from which we started on. Model is based on ordinary differential equations for species abundances. Main modifications from LL study concern (i) the limitation of the feeding kernel to only smaller species that we removed, what allows loops and cannibalism (ii) the functional response that is Holling-type II in our model, (iii) a simplified equation for the basal resource, and (iv) the introduction of two different feeding conversion coefficient, for the consumption of the basal resource or for the feeding on species. For convenience, we denote -→ ϕ i = x i , s i , d i ∈ R 3 the trait vector, thus the phenotype of species i. Following this notation, the feeding kernel of species i is denoted γ(x, -→ ϕ i ) where x is the body size of the prey. N i (t) denotes the biomass of species i at time t.

N 0 (t) is the biomass of the basal resource.

N(t) = N 0 (t), N 1 (t), • • • , N n (t) is the biomass vector of the community. Φ = -→ ϕ 1 , • • • , -→ ϕ n = (x 1 , s 1 , d 1 ), • • • , (x n , s n , d n )
is the trait vector of the community. The population dynamics reads :

                               . . . dN i dt = N i f P F P ( -→ ϕ i , Φ, N)
basal resource consumption

+ f C F C ( -→ ϕ i , Φ, N) species consumption -m( -→ ϕ i ) mortality -C(x i , Φ, N) competition -µ P (x i , Φ, N)
mortality from predation . . .

dN 0 dt = I -eN 0 -µ P (0, Φ, N)N 0
(1) Table 1 gives the detailed formulations of the different terms of these differential equations. The choices in the parameter values are based on the existing literature. Parameters values for the resource dynamics are based on LL model, as well as competition and mortality functions parameters. As in [START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF] we replaced the body size dependent efficiency function by the two fixed assimilation efficiency parameters, one for basal resource, f C , and one for species resource, f P , following [START_REF] Yodzis | Body Size and Consumer-Resource Dynamics[END_REF]. Parameters values are indicated in the caption of Figure 2. For the range of parameters values we considered, the joint dynamics of population sizes generally stabilizes, after a short transient dynamics, to a stable steady state.

Canonical equations. The invasion fitness can be inferred from Eq. ( 1) :

N(t) = N 0 (t) basal resource , N 1 (t), • • • , N n (t)
species Community traits: body size ; feeding kernel width ; prefered prey body size difference :

Φ = -→ ϕ 1 , • • • , -→ ϕ n = (x 1 , s 1 , d 1 ) species 1 , • • • , (x n , s n , d n ) species n
Gain from consumption on basal resource for species i :

F P ( -→ ϕ i , Φ, N) = γ(0, -→ ϕ i )N 0 1 + h( -→ ϕ i ) ∑ n k=1 γ(x k , -→ ϕ i )N k + γ(0, -→ ϕ i )N 0
Gain from predation on species for species i :

F C ( -→ ϕ i , Φ, N) = ∑ n j=1 γ(x j , -→ ϕ i )N j 1 + h( -→ ϕ i ) ∑ n k=1 γ(x k , -→ ϕ i )N k + γ(0, -→ ϕ i )N 0
Feeding kernel for species i on basal resource :

γ(0, -→ ϕ i ) = γ 0 e - (x i -d i ) 2 s 2 i
Feeding kernel for species i on species j :

γ(x j , -→ ϕ i ) = γ 0 e - (x i -x j -d i ) 2 s 2 i Handling time : h( -→ ϕ i ) = h 0 + h 1 s 2 i
Loss from predation for species i :

µ P (x i , Φ, N) = n ∑ j=1 γ(x i , -→ ϕ j )N j 1 + h( -→ ϕ j ) ∑ n k=1 γ(x k , -→ ϕ j )N k + γ(0, -→ ϕ j )N 0
Consumption rate for basal resource :

µ P (0, Φ, N) = n ∑ j=1 γ(0, -→ ϕ j )N j 1 + h( -→ ϕ j ) ∑ n k=1 γ(x k , -→ ϕ j )N k + γ(0, -→ ϕ j )N 0
Competition on species i :

C(x i , Φ, N) = n ∑ j=1 α 0 e - (x i -x j ) 2 σ 2 compet N j
Intrinsic mortality of species i :

m( -→ ϕ i ) = m 0 x -1 /4 i 1 1 -x i/x max + s 0 s i -1 /4 f ( -→ ϕ mut , Φ) = f P F P ( -→ ϕ mut , Φ, N(Φ)) + f C F C ( -→ ϕ mut , Φ, N(Φ)) -m( -→ ϕ mut ) -C(x mut , Φ, N(Φ)) -µ P (x mut , Φ, N(Φ)) (2)
where N(Φ) is the vector of equilibrium values of the populations biomasses, determined by the traits vector Φ.

From the invasion fitness, Adaptive Dynamics toolbox let us infer the canonical equation for each coevolving species and for each trait x, s and d. For species i :

dx i dτ = k x Ni (Φ) ∂ f ∂x mut ( -→ ϕ mut = -→ ϕ i , Φ) ds i dτ = k s Ni (Φ) ∂ f ∂s mut ( -→ ϕ mut = -→ ϕ i , Φ) dd i dτ = k d Ni (Φ) ∂ f ∂d mut ( -→ ϕ mut = -→ ϕ i , Φ) ( 3
)
where τ is the evolution time and k x , k s and k d are mutation-related constants. This system of 3n ODE drives the coevolution of the n interacting species in the community. Evolutionary dynamics usually stabilizes to an equilibrium of traits values, an evolutionary singularity. But evolution does not necessarily stops here. One of the power of Adaptive Dynamics is that it can predicts, again in the basis of the invasion fitness function, the two subsequent alternative outcomes : long term stabilization of traits value, or the occurrence of the evolutionary branching that leads to the emergence of two new species in the place of their unique parent species.

Trophic network emergence. We proceed with an initial single species (n = 1) for which we derive and integrate the canonical equations Eq. (3) for its three traits.

If the evolutionary dynamics stabilizes to an evolutionary singularity that allows evolutionary branching, we replace the former species with two new species with slightly different traits compared to their parent species. We rewrite Eq. ( 3) with now 2 species (6 equations) and integrate it again. We stop integration upon two possible events. (i) one of the species undergoes evolutionary branching: we replace it by two new species (ii) traits in the population co-evolve to a point where one of the population vanishes at equilibrium: we remove the corresponding species from the community. We rewrite Eq. ( 3) with the new number of species (n + 1 or n -1). We iterate the process until a chosen final evolutionary time is reached. A typical simulation is given by Figure 1 (panel II). This iterative approach of Adaptive Dynamics canonical equation has, to our knowledge, only been used but in a very recent study [START_REF] Hui | Modelling coevolution in ecological networks with adaptive dynamics[END_REF]. Compared to LL, this simulation method is much more efficient: simulations run faster and by construction, species are well separated from each other so that the final states of simulation can easily been analyzed. As already put out in [START_REF] Allhoff | Evolutionary food web model based on body masses gives realistic networks with permanent species turnover[END_REF], coevolution of the feeding kernel along with the body size allows the emergence of much more diverse trophic networks compared to models in which only body size evolves. On the other hand, as a downside, this additionnal freedom we inserted in the model released the individual parameters values (feeding kernel parameters d and s but also body size x ) from the intermediate range they were confined into in LL's model, so that we could not, as they did, neglect some of the constraints unveiled when parameters attain extreme values. We hence incorporated in the model additionnal costs due to high body size (high x) and to extreme specialization (very low s). Moreover, we modeled the handling time in the Holling-type II functional response as an increasing function of s so that a complete generalist consumers has a broader prey niche, but is less efficient in feeding than specialists. See table 1 for a precise formulation of these costs.

Analysis of evolving food webs in a spectrum of contrasted environmental conditions

Method. Trophic networks generated by evolution are highly dependent on the model parameters. Some of them are related to individuals or to interactions between individuals. Others are directly related to environment through the dynamics of the basal resource. We decided to focus our analysis on the global scale parameter representing the environmental richness, I (influx level of the basal resource, Eq. 1). This choice was first motivated by the fact that it allows to a certain extent to reproduce the contrasting environmental conditions of natural ecosystems. In addition, preliminary study shown that this parameter by itself allowed to generate a great diversity in simulated evolved food webs.

We explored some different values for the other parameters but then we fixed them.

For each value of I in a log scaled range of 1000 values from 0.1 to 350, we started from a single species, which traits (x 1 , s 1 , d 1 ) were randomly chosen, and we let our algorithm generate an evolved, diversified, food web. We let the simulations run for 4.10 6 time steps, what always allows the network to stabilize to some steady evolutionary final state. Then we analyzed each evolved trophic network for a series of descriptors:

• species diversity S (final number of species),

• connectance L /S(S+1), where L is the number of realized links,

• maximum trophic level TL max , computed as in (30), such that

TL i = 1 + S ∑ j=0 γ j,i TL j ,
• functional diversity indices for species attributes (body size x, generalism s, preferred prey body size xd) ; simplifying formula from (31), we computed these indices as the weighted variance of the attribute:

FDvar x = 1 N tot S ∑ j=1 N j (x j -x) 2
.

Figure 1 (panels I, II and III) synthetizes our method.

Results and Discussion. We compiled data from 213 natural food webs from the ECOWEB database [START_REF] Cohen | Ecoweb 1.1: Ecologists' cooperative web bank[END_REF], and determine the properties of these food webs with the purpose to compare them to the properties of our theoretical emergent food webs. Figure 2 (A-D) con-fronts our theoretical networks to these natural food webs for the relationships between their specific diversity and their connectance, maximum trophic level, percentage of "top" species (consumed by no other species) or percentage of "bottom" species (that feed only on basal resource). The general trends we found were similar. Besides, although we only varied a single control parameter, the relationships between these emerging properties of evolved food webs are not curvilinear, the dots clouds really spreading over a 2D region. Even if natural data usually better fill the range of variation of some properties, we were able to reproduce a significant part of the diversity of trophic structures observed in the wild. This could suggest that the variability observed in the natural trophic networks might be explained by a limited number of factors.

We have not done more precise statistical tests here because it does not make much sense for several reasons. First of all, mere data fitting can be a poor test of model performance in many cases in ecology [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF][START_REF] Mcgill | Rebuilding community ecology from functional traits[END_REF] In ad-

α 0 = 0.5, σ compet = 0.4, h 0 = 0.5, h 1 = 0.05, γ 0 = 1, m 0 = 0.1, x max = 15, s 0 = 0.1, f p = 0.85, f c = 0.45.
dition, our species and their traits are not readily reliable with this set of empirical data as they do not consider body size and often pool several species in one node. Moreover, in empirical food webs, data are all the same compiled over a wide range of environmental situations, making it difficult to directly compare the variations of characterizing properties with a single environmental determinant as the environmental richness we focused on. Empirical studies that examine the links between such large-scale parameters and community structures do exist, though [START_REF] Kortsch | Food-web structure varies along environmental gradients in a highlatitude marine ecosystem[END_REF]. Figure 3 (A-C) shows in more details the variations in selected properties of emerging food webs against the gradient of environmental richness we considered. We chose to emphasize functional diversity indices for the adapting traits of species in the networks. The first striking result is that these indices of functional diversity, depending on which trait we measured them, shape differently. Functional diversity index on maturity body size x increases all along with environmental richness. Emerging specific diversity also follows the same trend, that is a monotonous increase with environmental richness (result not shown). Functional diversity index on generalism level s and on prefered prey body size xd as for them are both maxized for intermediate values of the control parameter I. But the values of I that maximizes the functionnal diversity of s are low to intermediate whereas those maximizing the functionnal diversity for xd are distinctly higher. Again, even if they do not appear on the figure, we additionnally tested connectance, which curves follows the shape of the functionnal diversity of s whereas the maximum torphic level follows the shape of the functionnal diversity of xd. The interesting point here is that these results make it possible to formulate hypotheses as to the key role of our control parameter, and only it, in the resulting emergent structure of food webs. In particular, three highlighted properties can not be maximized at the same time, because they are maximized when trophic networks evolved in three different environmental richness values. Highest specific diversity does not necessarily implies highest trophic complexities. This last point was previously highlighted in [START_REF] Brännström | Emergence and maintenance of biodiversity in an evolutionary food-web model[END_REF]. Besides, previous studies have also investigated the role of environmental attributes on food web structure, and especially on food chain length. Three hypothesis have been formulated [START_REF] Post | Ecosystem size determines foodchain length in lakes[END_REF], namely the productivity hypothesis, the ecosystem size hypothesis, and the combined productive-space hypothesis. While this pioneering study shown a strong inclination for the ecosystem size hypothesis, other studies have shown results in favor of the other two. Our result on the complex structuring role of environmental productivity can give support to this productivity hypothesis. An advantage of our study is that we could isolate the effects of this structuring factor and test variations on a wide range of theoretical values. To illustrate this, we can compare our results with a recent study that looks at structural variations in marine communities at high latitudes along an environmental productivity gradient [START_REF] Kortsch | Food-web structure varies along environmental gradients in a highlatitude marine ecosystem[END_REF]. Like them, we find the lowest connectance values for low environmental productivity, and higher connectance values, or higher trophic levels for high productivity values. However, in our case, the trophic complexity (connectance or functionnal diversity) of evolved food webs, after a first increase to some maximum value then begins to decrease again for high environmental productivity. This underline the fact that relationship between global scale parameters and network structure are not always straightforward and monotonuous and could give cues of why whichever simple hypothesis may be verified in some situations but defaulted in others.

Disturbance experiment

Method. There are numerous ways to understand and quantify the stability of an ecosystem [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF]. A large part of them consist in quantifying the resistance to perturbations. In this disturbance experiment, we exposed evolved food webs to new environmental conditions, different from the acustomary level of environmental richness during their evolution (Figure 1, all panels). Effects of perturbations are assessed on a short time scale : once a food web has reached its evolutionary steady state, we abruptly change parameter I in the demographic system [START_REF] May | Unanswered questions in ecology[END_REF]. It may only result in a change in the population abondances that stabilize, on the demographic timescale, to some modified equilibrium values. But it often happens that in this new equilibrium, one or several of the populations vanish, issueing species losses in the trophic network. We can therefore estimate how evolved food webs respond for various, new conditions, in terms of species diversity change (assessing species loss). But all other properties of food webs are likely to be impacted as well. We chose to focus on two properties. The maximum trophic level property, as in the previous section, and a new "information characteristic": entropy. Entropy, or "information diversity", is derived from Shannon's information index [START_REF] Shannon | The mathematical theory of information[END_REF], that quantifies here the diversity of trophic exchanges in a network [START_REF] Macarthur | Fluctuations of Animal Populations and a Measure of Community Stability[END_REF][START_REF] Rutledge | Ecological stability: An information theory viewpoint[END_REF][START_REF] Robert E Ulanowicz | Biodiversity, functional redundancy and system stability: Subtle connections[END_REF]. It is computed as:

H = - n ∑ i n+1 ∑ j T ij T .. log T ij T .. , (4)
where T ij denotes the magnitude of flow from node j to node i. Each food webs, evolved in a richness level I = I evol , are confronted to a range of 1000 new richness levels I perturb from 0.1 to 350. We hence get a 1000 × 1000 treatment. For each configuration, we determine the relative species loss, that is the number of perturbationmediated extinctions, related to the pre-disturbance number of species. On contrasted situations, we assess the robustness of food webs, defined as the maximum intensity of disturbance they can undergo before they lose more than half of their species. We also calculate the gain in maximum trophic level, which can be positive or negative. It is obviously difficult to determine a posteriori the environmental conditions an ecosystem has evolved in. This first approach allows us to determine how ecosystems react in face of changes in their environmental richness level, but we found the necessity to link an observable property of contemporary evolved food to their fragility or robustness. In a second approach, we chose to relate the connectance of evolved food webs to their relative species loss when confronted to a reduction factor in their richness level : I perturb /I evol ranging from 0.2 to 1. As the connectance is an emergent property and not a control parameter, we used a smoothing gaussian kernel to merge data from irregularly distributed connectance values. The local variance of data we averaged with the smoothing kernel varies significantly and we quantified it precisely. Lastly, we interested in the comparison of a trophic network that evolved in a given level of richness I evo to other networks that evolved in other environmental conditions but were then subjected to this richness level as a perturbation. We compute the comparative species loss, corresponding to the difference between the number of species of a trophic network perturbed by a new richness level I and the number of species of a food web evolved from the beginning in this richness level I. Similarly, we computed the comparative entropy loss.

Robustness of evolved food webs. Figure (4,A) shows the relative species loss for all 1000x1000 treatment configurations. It highlights three main outcomes. Overall, food webs appear to be quite robust to moderate environmental changes : no catastrophic collapses, relatively low relative species losses. Second, evolved food webs appear to be much more sensitive, in terms of species loss, to environmental depletion than to enrichment. Enrichment leads to a maximum of 10 -15% species loss and often does not have any impact on specific diversity. Lastly, we found that food webs that have evolved in poor environments are less robust to environmental impoverishment than webs that have evolved in a rich environment. For a similar reduction factor in the environmental richness, food webs evolved in a rich milieu, even in terms of relative species loss, suffer much less of the perturbation than food webs evolved in poor conditions. The red arrows show the disturbances necessary to reduce the number of species by half, for two contrasted I evo . It happens for a factor 10 of reduction of then richness level for the left arrow whereas it needs a factor 200 for the right arrow. When considering the maximum trophic level (Fig. 4,B), it appears, contrary to the results on species loss, that the evolutionary history parameter I evo is not so significant any more, but that the maximum trophic level always tends to be maximized for some intermediate value. Indeed, trophic networks evolved in poor conditions increase their maximum trophic level when they are immersed in richer conditions, whereas for those evolved in rich environment, it increase when environment is impoverished. Now relating the emerged connectance of evolved food webs back to relative species loss (Fig. 4, C), we focused on environment depletion. We observed that evolved trophic networks characterized by low connectance values appear to be more robust to an impoverishment of their environmental condition. Even for a division by 5 of the environmental richnes, they only lose 20% of their species. By contrast, the most connected of our evolved food webs lose about half of their species for only a division by two of the environmental richness. Yet, we assessed that the variance of the averaged species loss becomes important, espcially for low connected food webs experiencing strong per-turbations (bottom left of Fig. 4,D). This results from the fact that a variety of evolved food webs, that experienced very different I evo during their evolution, end up with convergent connectance values. Even if these low connectance food webs are on average much robust to perturbations, some of them, from our panel, evolved in poor environmental conditions and are highly impacted by further impoverishment. Both topological and dynamic approaches were used to study community responses to species loss [START_REF] Dunne | Network structure and biodiversity loss in food webs: Robustness increases with connectance[END_REF][START_REF] Solé | Complexity and fragility in ecological networks[END_REF][START_REF] Eklöf | Species loss and secondary extinctions in simple and complex model communities[END_REF]. Dynamic approaches were used, taking into account some of the weaknesses of topological approaches (e.g. misleading picture of potential secondary extinctions). However, most of these are based on generalized LV models, ignoring the evolutionary history of the communities generated. Explicitely taking into account this evolutionary history, we get in our study a data set arguably closer to field observations. The second type of divergence concerns the type of disturbance. These studies were carried out through the analysis of secondary extinctions following the initial loss of one or more species. Primary species losses are discussed as being the result of several disturbance factors, that however usually stay out of the scope of the studies. Even if the results may vary depending on how the species are primarily removed (e.g. primary or secondary consumers), one of the major findings of these studies is that food webs characterized by high connectance values appear to be the most robust.

Here we have a different approach. When primary extinctions are proximal causes for secondary extinctions, explicitly taking into account changes in environmental conditions is like considering distal causes for these extinctions. And hence, we do not decide a priori which species potentially disappear but we estimate it a posteriori from the new environmental parameter value along with the population dynamics model. And therefore, our apparently opposite results, because they were obtained on a selected data set and with a different methodology, actually do not contradict the current paradigm around the connectancerobustness debate. Instead, we provide additional elements of understanding and underline the precaution of attributing a high robustness to highly connected food webs if the environment is to undergo depletion disturbances. We also highlight that food webs with similar level of connectance may respond very differently to perturbations so that the average response may not always be very informative. This final result is obviously very sensible on the data set used and we recommend that precautions be commonly taken. Another advantage of the framework used here is that it allows us to look at the impacts of an explicit disturbance, here the change in environmental richness, on other aspects of biodiversity than specific diversity. In particular, studying trophic level loss/gain can be allowed because we have access to the demography of the different species, which is not the case with topological approaches for example. In our setting, perturbations do not change trophic connections between species. The results where maximum trophic level increases due to perturbations could only be explained because demography was rebalanced. For example, for cases where maximum trophic level increases with richness reduction, biomass of low level species, that exclusively consumes basal resource, decreases, thus representing a reduced proportion of the diet of higher level consumers, whose trophic level therefore increases.

Adaptation of evolved food webs. Results on maximum trophic level (Fig. 4,B) shows that evolutionary history do not always shape food webs in a way that their functionning, at least for some characteristics, be optimized for the environmental conditions that have accompanied its evolutionary emergence. Indeed, all food webs increase their maximum trophic level when they are immersed in some optimal intermediate level of richness I opt , regardless of their evolutionary conditions. Undoubtedly, a food web is not an adaptive unit that directly undergoes natural selection. Still, we wanted to push this analogy and assess whether trophic networks could be designed by co-evolution so that they would particularly fit their environment for some other property. With that goal, we observed our setting as a common garden experiment: we compared food webs "transplanted" in a given environment, one of them evolved in the same conditions, the second one in an other condition, richer or poorer. Figure 4 (E) shows the comparative species loss of the "exogenous" food web compared to the "native" food web. The striking result is that, concerning specific diversity, food webs that emerged in poor or intermediate level of richness are not especially best "adapted" to their environment. Indeed, most of the trophic networks that evolved in a richer milieu, and then gained a greater specific diversity, keep a higher number of species when they undergo even severe environment impoverishment, compared to food webs evolved from the beginning in these poor environments (red area under the diagonal). Figures (4,E,F) allows us to estimate whether the evolutionary history of a food web leads it to be more adapted to its evolutionary environment than any other food web that undergoes these conditions after its emergence. With regard to specific diversity, it appears that food webs that have evolved in rich environments and are then confronted to an impoverished environment have more species than food webs that have adapted to this poor environment during their evolution. It therefore appears that evolved trophic networks do not seem at first sight to be the most suitable for their evolving environment, at least not in view of the specific diversity present. For connectance, not presented results show a symmetrical situation: networks that have evolved in poor environments, which then subsequently enriched, have (still) higher connectance than networks that have evolved in these rich environments. It appears that the conclusions of which is best fitted to a given environment thus may vary according to the properties studied. However, when looking at properties derived from information theory, we found that evolution could optimize food webs for some characteristics. Indeed, the properties related to the organization and diversity of flows within evolved trophic networks show this trend. For most evolved food webs, any disruption of environmental richness from accustomary conditions induces values of these properties that are lower than those that characterize food webs that have evolved in these environments (Fig. 4,F). This means that the structure of food webs and the flows within them seem optimized by evolving processes. The fact that ecosystems may develop or evolve by maximizing some of their properties is a long, contentious debate in ecology [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF]. The approach discussed here allows us to question the level of adaptation of a food web to its evolutionary environment by comparing it to other food webs subjected to this environment without an adaptation period. This provides a complementary point of view to the extent to which evolution leads an entity, here, in this case, a food web, to develop optimal properties. Considering food webs as super-organisms is a concept widely discussed in the evolutionary and ecological literature [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF][START_REF] Sloan | Reviving the superorganism[END_REF]. Loreau [START_REF] Loreau | From populations to ecosystems: theoretical foundations for a new ecological synthesis[END_REF] argues that while near optimization of ecosystems properties is conceivable (with constraints on strengths and directions of individual and group selections), there is no reasons to believe that it will be generally achieved, and that there is no objective criterion to forecast it. Here we present contrasted results. Contrary to what usually show cross-transplantation experiments on real organisms, for food webs, some of them are better fitted to some environments than other food webs "adapted" to them (results on number of species, connectance). Besides, ecosystem properties related to diversity and organization of flows appear to be the best candidates for a generalized near-optimal behavior.

Conclusions

In this paper, we first focus on the influence of a global scale parameter, environmental productivity, on the structure of food webs. Using a community evolution model, we look at the extent to which variations in this control parameter alone can lead to the emergence of food webs that reflect the diversity of structures observed in the wild. We show that a considerable part of this diversity can be reproduced in this way, underlining the importance of environmental productivity as a determinant of community structure. The second particularity of this study is to study the shortterm ecological responses of the food webs produced in this way. We focus on an explicit disturbance mode, which aims to modify the control parameter once the evolution step is completed, and thus quantify a posteriori the losses (or gains) of different network properties. Food webs characterized by high connectivity values appear to be the least robust to a decline in their environmental condition, given the resulting percentages of species loss. The induction of intermediate environmental productivity values positively (or neutrally) affects all food webs, with respect to the maximum trophic level encountered. This disturbance experiment finally allows us to question the potential status of food webs as "super-organisms", i.e. their state of adaptation to their evolving environment, by comparing and contrasting other food webs subject to this environment. Some properties, linked to the diversity of flows within networks, appear to be the best candidates for possible properties that would thus be maximized by evolutionary processes. We believe that this work represents a first step in the unexplored avenue of exploiting community evolution models to analyze responses to disturbances of such winnowed survivors.

Fig. 1 .

 1 Fig. 1. Explanatory scheme of the disturbance experiment. The first step consists in setting up the model for different environmental richness values I (panel I). We show here three example values : poor, intermediate and rich environmental conditions. The second stage (panel II) consists in simulating the eco-evolutionary model under these environmental conditions. Trophic networks emerge from succesive diversifications. Evolved food webs are analysed (panel III) for different emerging properties. Then (panel IV), they are confronted to changes in environmental richness, what impact populations on a short time scale, potentially inducing species extinctions (black circles).

Fig. 2 .

 2 Fig. 2.Emergent properties of evolved food webs (blue) compared to those of empirical food webs (red). Bottom corresponds to the percentage of species feeding on the basal resource. Top corresponds to the percentage of species that are not consumed. Parameters, expect I, are fixed and take the following values : α 0 = 0.5, σ compet = 0.4, h 0 = 0.5, h 1 = 0.05, γ 0 = 1, m 0 = 0.1, x max = 15, s 0 = 0.1, f p = 0.85, f c = 0.45.

Fig. 3 .

 3 Fig. 3. Functional diversity indices against the evolutionary environmental richness I. Functional diversity is computed for (A) body size, (B) degree of generalism, and (C) optimal prey body size.

Fig. 4 .

 4 Fig. 4. Short-term, ecological responses of food webs when they are subjected to new environmental conditions different from their evolutionary environment. (A,B) present respectively the relative species losses and the gains/losses of maximum trophic level observed. (C) shows the relative species losses when the productivity of the basal food web resource is depleted by a certain percentage. Food webs are sorted by connectance, and values are smoothed by a Gaussian kernel. (D) shows the local standard deviations of relative species losses. (E,F) show, respectively, the comparative losses of specific diversity and entropy (diversity of flows) of evolved trophic networks. Comparative loss refers to the difference between the new property value of a food web for a given, new, environmental richness, and the of the same property for the "native" food web that has evolved in this later environmental richness.

Table 1 .

 1 formulasPopulations size :

Girardot | Does evolution design robust food webs ?

Acknowledgments

We thank Clement Aldebert for helpful discussions and comments. This article is part of the CIGOEF French ANR. The PhD scholarship of BG was funded by the French Ministry for Education and Research.

Authors contribution

B.G. and M.G. designed the model, performed the study and analysed the results. All authors contributed to the writing of the manuscript. The authors declare no conflict of interest.