A Gutierrez

E Marcault

C Alonso

J-P Laur

D Tremouilles

Parseval's Theorem for the Inductor Analysis in High-Frequency Power Converters

Keywords: Power Inductor, Quality Factor, Self Resonance Frequency, Boost Converter, GaN-HEMT

Advances in GaN-HEMT devices have encouraged the development of power converters at megahertz-level. However, inductors for high-frequency power converters require innovative approaches to overcome challenging constraints in both power and frequency. In this context, this paper presents an extension of the Parseval's theorem for the inductor analysis based on the energy conservation principle in the time and the frequency domains. Our approach aims to provide insights about the inductor power losses using high-frequency parameters and the inductor current harmonics. The proposed approach disaggregates the power losses in the frequency-domain for the inductor power signals in the time-domain. Main findings from the presented methodology provide useful selection criteria for inductors in power converters given parameters of quality factor (Q) and Self Resonance Frequency (SRF ). Simulation results show the impact of the inductor behavior on the switching losses. An experimental setup validates the proposed approach. A high-frequency boost converter is presented as a study case.

Introduction

Traditional analysis of power converters requires design tools intended to deal with operation frequencies around hundreds of kilohertz [START_REF] Ning | Review of si igbt and sic mosfet based on hybrid switch[END_REF]. However, outstanding advances in Gallium Nitride -High Electron Mobility Transistors (GaN-HEMTs) enlarge the operation condition of power converters to the range of megahertz [START_REF] Kashiwagi | Hysteretic control embedded boost converter operating at 25-mhz switching[END_REF] [START_REF] Hariya | Five-megahertz pwm-controlled current-mode resonant dc-dc step-down converter using gan-hemts[END_REF]. Indeed, this extended bandwidth requires innovative design methodologies to achieve a trade-off between high-frequency and high-power conditions [START_REF] Rodríguez | High-frequency pwm buck converters using gan-on-sic hemts[END_REF]. Additionally, these new generation of power converters should incorporate power inductors with both high current handling and high operation frequency [START_REF] Hilal | Optimal inductor design and material selection for high power density inverters used in aircraft applications[END_REF] [START_REF] Liang | 3-d-printed air-core inductors for high-frequency power converters[END_REF].

As result, high-frequency parameters become more and more important for the analysis and selection of inductors during the design process of power converters [START_REF] Jaimes | Characterization of high-q inductors up to its self-resonance frequency for wireless power transfer applications[END_REF] [START_REF] Corrêa | Experiments with a compact wireless power transfer system using strongly coupled magnetic resonance and metamaterials[END_REF]. Usually, applications of inductors in high-frequency take into account the manufacturer parameters of Self Resonance Frequency (SRF ) and quality factor (Q). These SRF and Q parameters allow describing the inductor operation bandwidth and the associated power losses. However, these inductor parameters are few analyzed in the conventional design of power converters at kilohertz-level. Conversely, operation at megahertz-level should include the SRF and Q parameters in the design stage.

Additionally, the switching improvements provided by GaN-HEMTs can extend the operation frequency of power inductors to levels usually dedicated to radio-frequency applications. However, manufacturing constrains associated to ferromagnetic cores limit the scope of power inductors in highfrequency applications [START_REF] Yao | A novel windingcoupled buck converter for high-frequency, high-step-down dc-dc conversion[END_REF]. Therefore, manufacturing features and powerfrequency behavior require a suitable trade-off to improve the inductor performance in high-frequency power converters [START_REF] Barrios | Optimal dc gapped inductor design including high-frequency effects[END_REF] [START_REF] Barrera-Cardenas | Optimal design of aircore inductor for medium/high power dc-dc converters[END_REF].

In this context, we proposed an innovative approach using the Parseval's theorem as a means to associate the power losses in the time-domain and the inductor behavior in the frequency-domain. The developed framework aims to be complementary to conventional and well known methodologies for the inductor selection and design of power converters. This approach provides insights about the impact of the quality factor (Q) and the Self Resonance Frequency (SRF ) on the inductor current harmonics and the power losses.

Simulation results suggest suitable criteria for a trade-off between the quality factor (Q) and the power losses given feasible manufacturing features. Furthermore, this paper proposes an inductor model incorporating the Q and SRF parameters into a circuital model useful for the analysis of power converters.

This paper is organized as follows. Section 2 describes the power and fre-quency behavior of power inductors considering high-frequency parameters. Section 3 presents the theoretical approach based on the Parseval's theorem to study the relation between power and frequency. Section 4 describes the proposed inductor model including SRF and Q parameters. Section 5 explains in detail the experimental setup to validate the proposed approach.

Finally, an experimental high-frequency boost converter is implemented.

Power inductor behavior in high-frequency

This section describes the inductor behavior in a high-frequency power converter taking into account the quality factor (Q) and the Self Resonance Frequency (SRF ). A high-frequency boost converter is presented as an illustrative example.

Fig. 1 shows the designed boost converter as a study case using a commercial reference of a power inductor. The switching frequency is set to 30MHz to take advantage of the switching characteristics of GaN-HEMTs. This boost converter increases the voltage from 200V to 400V with output power of 400W and load of 400Ω. In this design, the inductor is set to 8.2µH and the output capacitance is set to 220pF. Furthermore, the studied converter focuses on the inductor performance considering ideal the other components to avoid their influence in the developed analysis. The power inductor under study accomplish the requirements of inductance and current in the range of KHz. However, it has a low performance in the range of MHz as shown in Fig. 2. The proposed analysis intentionally begins with this unsuitable inductor to assess and understand the influence of the high-frequency parameters in the inductor performance. This analysis is the first step to define selection criteria for high-frequency power inductors with a suitable trade-off between performance and feasibility. Fig. 3 shows the simulation results comparing the inductor current in the case of an ideal inductor in series with a low resistance and the manufacturer model for the studied inductor. In addition, Fig. 4 illustrates the current harmonic spectrum of the inductor under analysis. Results from Fig. 2 to Fig. 4 show that the unsuitable SRF causes distortion in the inductor current given the harmonics in the capacitive region beyond of the SRF of 22MHz. In addition, the low Q parameter causes important power losses in comparison with the ideal inductor in series with a low resistance. These power losses in the time-domain are associated to the interactions between the inductor harmonics and the inductor behavior in high-frequency. Next section will discuss these interactions using the Parseval's theorem, which describes the energy conservation in the frequency-domain and the time-domain.

Parseval's theorem approach

As described in the previous study case, the wrong selection of the power inductor leads to low signal quality and higher losses in high-frequency power converters. As a result, the following theoretical approach provides insights about the relation between power and frequency of inductors suitable for high-frequency power converters.

Considering the approach described in [START_REF] Kelkar | An extension of parseval's theorem and its use in calculating transient energy in the frequency domain[END_REF] and defining the inductor current i(t) and the inductor voltage v(t), the electrical energy U L for an inductor in the time-domain is given by eq.( 1).

U L = ∞ -∞ p(t) dt = ∞ -∞ i(t)v(t) dt (1)
In the frequency-domain, the convolution theorem applied to the power expression i(t)v(t) of eq.( 1) is given by eq.( 2),

F {i(t)v(t)} = I(f ) * V (f ) (2)
where I(f ), V (f ) are the complex valued Fourier transforms. F denotes Fourier transform, and * denotes convolution. By definitions of Fourier transform F and convolution, the eq.( 2) becomes eq.( 3),

∞ -∞ i(t)v(t) e -j2πσ dt = ∞ -∞ I(f )V (σ -f ) df (3)
the evaluation of the Fourier transform at the origin (σ = 0) equals the integrals over all their domains [REF REF]. Thus, the electrical energy U L from eq.( 1) can be expressed by eq.( 4), where V (f ) is the complex conju-

gated of V (f ), ∞ -∞ i(t)v(t) dt = ∞ -∞ I(f )V (f ) df (4) therefore, ∞ -∞ i(t)v(t) dt = ∞ -∞ I(f )[I(f )Z(f )] df = ∞ -∞ I(f ) 2 Z(f ) df (5)
given that Z(f ) = R f + jX L is the complex impedance of the inductor and X L is an odd function [START_REF] Kelkar | An extension of parseval's theorem and its use in calculating transient energy in the frequency domain[END_REF],

∞ -∞ i(t)v(t) dt = ∞ -∞ |I(f )| 2 (R f -jX L ) df = ∞ -∞ |I(f )| 2 R f df (6)
the result in eq.( 6) illustrates the energy conservation between the time and the frequency domains. Therefore, eq.( 6) can be seen as an extension of the Parseval's theorem [START_REF] Kelkar | An extension of parseval's theorem and its use in calculating transient energy in the frequency domain[END_REF]. The physical interpretation of eq.( 6) is that the total energy U L of the inductor can be calculated by integrating power over time or by the spectral power across frequency known the current and the associated losses. The Parseval's identity define the relation between the average power of a signal h(t) and their Fourier coefficients as,

1 T T /2 -T /2 |h(t)| 2 dt = ∞ n=-∞ |C n | 2 (7)
therefore, the average inductor power can be expressed from eq.( 6) by, 8) is applied to the designed boost converter by means of the associated inductor losses R f and the current harmonics (see Fig. 2 and Fig. 4). The results for the average power in the time and frequency domains are summarizes in Fig. 6 and Fig. 7. These results agree with the expected energy conservation criterion.

P AV G = 1 T T 0 i(t)v(t) dt = N n=0 I 2 n R f (8)
Fig. 7 depicts the distribution of power losses in the frequency domain.

Results in Fig. 7 show that harmonics higher than eight times the switching frequency F sw have a negligible impact on the power losses. In addition, this figure allows highlighting the contribution to the power losses of the current at the switching frequency. In this case, the higher losses are given at the switching frequency F sw despite of the very high ratio between the DC current and the current at F sw . The relative high value of R f explains this power losses at the switching frequency F sw . Therefore, it is fundamental to increase the quality factor Q to decrease the R f losses at the switching frequency F sw in order to improve the global inductor efficiency.

Fig. 8 shows the relation of the power losses in the time domain against the quality factor at the switching frequency Q(F sw ),

Q(F sw ) = X L (F sw ) R f (F sw ) ( 9
)
The plot is calculated from eq.( 8) by assuming that the highest contribution to the power losses comes from the fundamental frequency at the switching frequency. To plot this figure, the SRF is assumed to be eight times the switching frequency F sw and the Q factor is evaluated at 30MHz for a fixed current spectrum. Results in Fig. 8 allow concluding that increasing considerably the Q factor has low impact in the power losses since the reduction in power losses becomes negligible. Therefore, it is necessary an approach to the suitable selection of the Q factor. In this context, the next section will propose a circuital model to include the SRF and Q parameters in the converter design process to assess their impact in the converter performance and to further inductor selection or manufacturing.

Proposed inductor model including SRF and Q

Currently, development of power electronics requires power inductors able to operate in high-frequency with high current capabilities. Therefore, the manufacturing specifications of SRF and Q should be useful to both power converter designers and inductor manufactures. As a result, this section provides a framework to integrate these parameters in the design process of power converters. For the inductor model of Fig. 9, the inductor impedance Z is given by the eq.( 10) at the switching frequency F sw and ω sw = 2πF sw ,

Z = 1 1 R p + 1 -ω 2 sw LC p + jω sw C p R s R s + jω sw L (10)
considering the real and imaginary parts of Z and solving for Q sw , where Q sw is defined as the quality factor at the switching frequency Q(F sw ), we have,

Q sw = Im(Z) Re(Z) = -R p (ω sw C p R 2 s -ω sw L + ω 3 sw C p L 2 ) R 2 s + ω 2 sw L 2 + R p R s (11)
solving for R p ,

R p = Q sw (R 2 s + ω 2 sw L 2 ) -ω 3 sw C p L 2 + ω sw L -ω sw C p R 2 s -Q sw R s (12)
Additionally, the SF R is defined by the inductor resonance frequency.

Thus, the C p capacitance is given by,

C p = 1 (2π) 2 (SRF ) 2 L (13)
Eq.( 12) and eq.( 13) allow including in the circuital inductor model the high-frequency parameters of Self Resonance Frequency (SRF ) and the quality factor at the switching frequency Q sw . Therefore, these expressions are useful to simulate and evaluate the impact of the SRF and Q parameters on the power converter performance.

As illustrative example, the boost converter of Fig. 9 is simulated for F sw =30MHz, L=8.2µH, and R s =0.2Ω. R p and C p are calculated from eq.( 12)

and eq.( 13). Q sw =100 in the analyzed cases. The simulation results of Fig. 10 show the impact on the drain current of the GaN-HEMT transistor when SRF is evaluated for 30MHz and 250MHz. These results depicts the increasing in the drain current by around 50% which lead to an increase in the switching losses by around 20%. This phenomenon is mainly generated by the increase of the parasitic capacitance C p when the SRF parameter is lower. The proposed analysis for power inductors has been described through this document. The simulated results have shown a suitable trade-off between power and frequency performance. Next section will introduce the experimental setup to validate the proposed approach.

Experimental results

This section describes the performed tests. First, the power inductors are measured and the circuital models are validated. Then, an experimental setup allows verifying the described power equivalence between the time and frequency domains. Finally, an experimental boost converter provides results about the performance of a conventional power inductor in a relative highfrequency.

Inductor circuital model and characterization

This test employs an impedance analyzer Agilent 4294A to measure the parameters of several power inductors. The impedance analyzer sweeps the frequency from 40Hz to 110MHz and it measures Z -θ, R -X and L -Q.

Table 1 summarizes the measured and calculated parameters. R p and C p calculated from eq.( 12) to eq.( 13) selecting the Q factor at 10MHz in all cases. Fig. 14 depicts the simulated and measured R f (the real part of the inductor impedance) for the inductor SRP5015TA. The MAPE (Mean Absolute Percentage Error) evaluates the model accuracy using eq.( 14) where M K is the measured value and S k is the simulated value. As listed in Table 1, the MAPE shows a partial agreement between the experimental data and the circuital model.

M AP E = 1 n n k=1 M k -S k M k (14)
However, Fig. 14 allows concluding that values for frequencies farther to the selected Q at 10MHz have less agreement than values around the selected Q. Therefore, the circuital model can slightly lose accuracy in a wide range of frequencies. The test is carried out at 10MHz and 30MHz. The aim of this test is to validate the power equivalence using the proposed approach. and case (b) for 30MHz for the power inductor SRP5015A. In the case (a), the fundamental frequency is lower than the SRF . Therefore, the inductor is able to store energy as a magnetic field with relative low power loss. In contrast, the case (b) has a fundamental frequency higher than the SRF . As result, the inductor behaves as a capacitor distorting the current signal and increasing the power loss.

Fig. 16 shows the FFT for the inductor current of case (a) in Fig. 15a.

In Table 2, the contribution of each harmonic is calculated from the measured R f (see Fig. 14) and the FFT (see Fig. 16) for the fundamental and harmonic frequencies using expression eq.( 15). Results from Fig. 15a and Table 2 agree with the expected correlation between the power in the time and frequency domains. Table 3 summaries the power results for the measured inductors. These results confirm the duality between the inductor power in the time-domain and the frequency-domain. However, differences between theoretical and experimental results are mainly caused by the shifting of the SRF given the parasitic capacitance and inductance of the current and voltage probes.

P AV G = 1 T T 0 i(t)v(t) dt = N n=0 I 2 n R f (15)

Boost converter at 1MHz

This test implements a boost converter at 1MHz. The boost converter specification are V in =30V, V out =60V, P out =40W, and inductor current I L =1.5A considering an efficiency η=0.9. The design uses and inductor SRR1210-8R2

given its favorable quality factor Q(1MHz)=30 and SRF =22MHz. The selected inductor SRR1210-8R2 has a series resistance R s =17mΩ and a saturation current of 7.5A. The switching frequency is set to 1MHz to ensure the eight harmonic lower than the SRF . A GaN-HEMT is used as a switching device. The achieved results harmonize with the theoretical framework.

However, the slight deviation between the temporal and frequency responses can be overcame by improving the experimental setup decreasing the parasitic elements. Fig. 17 shows the inductor waveforms. The inductor average power in the time-domain is P avg =492mW. Fig. 18 depicts the harmonic spectrum for the inductor current. Table 4 lists the power contribution of each harmonic to the total power. The total power in the frequency-domain is P avg =447mW.

Results from Fig. 17 and Table 4 validate the equivalence between the power in the time and frequency domains. In additions, the results allow identifying that the behavior at the fundamental frequency is the main cause of the inductor power loss. This section shows the results of implementing a high-frequency boost converter. A power inductor was designed for this experimental setup. The designed inductor has specifications of L=8µH, I sat =3.5A, Q(30M Hz)=25

and SRF =200MHz. The switching device is a GaN-HEMT.

The boost converter is tested first at V out =60V -P out =10W and after at

V out =400V -P out =400W. The test is carried out first in relative low power to use available current probe (Tektronix CT2) because current probes for high-current, high-voltage and high-frequency are currently under development. This represents a research challenge. In the second case, the inductor current is not measure for previous reasons. However, the global power converter results shows the suitable performance of the designed inductor.

In the first test, the specifications of the boost converter are V in =30V,

V out =60V, and P out =10W at 30MHz. Fig. 19 shows the waveforms for the power inductor. In this case, the average power was P avg =780mW in the time-domain and P avg =700mW in the frequency-domain.

Conclusions

The reported methodology has associated the analysis in the time and frequency domains for inductors in power converters using an extension of the Parseval's theorem. The proposed approach is complementary to conventional methodologies for the design of power converters. This analysis methodology allowed determining suitable criteria for the selection and simulation of inductors according to expected power losses. The proposed approach allowed disaggregating the power losses in the frequency-domain for complex inductor power signals in the time-domain. The study of the quality factor Q allowed concluding that increasing considerably Q has low impact on the power losses since the reduction in power losses becomes negligible.

Therefore, a moderate Q factor can be selected to achieve a trade-off between inductor performance and manufacturing feasibility. Additionally, the Self Resonance Frequency SRF around eight times the switching frequency is a suitable criterion to avoid inductor current distortion. The proposed model in this work included frequency parameters of inductors in the design process of power converters. The experimental results have validated the proposed approach. However, the experimental setup should decrease the parasitic inductance and capacitance to minimize the measurement disturbances at high-frequency.

Figure 1 :

 1 Figure 1: High-frequency boost converter with actual inductor and ideal associated components.

Fig. 2

 2 Fig. 2 depicts the frequency behavior of the associated inductor losses R f and the quality factor Q f for the power inductor of 8.2µH used in the designed boost converter. Additionally, Fig. 2 shows the SRF of the studied power inductor. The inductive characteristics prevail in the region of

Figure 2 :

 2 Figure 2: Frequency behavior of R f and Q f for the power inductor used in the designed boost converter.

Figure 3 :

 3 Figure 3: Inductor current. Color nomenclature: blue -ideal inductor in series with a low resistance, red -inductor under study.

Figure 4 :

 4 Figure 4: Current harmonic spectrum for the studied inductor.

Figure 5 :

 5 Figure 5: Extension of Parseval's theorem for the inductor analysis in power converters.

Figure 6 :

 6 Figure 6: Inductor power in the time-domain.

Figure 7 :

 7 Figure 7: Inductor power in the frequency-domain.

Figure 8 :

 8 Figure 8: Inductor losses and quality factor Q sw at the switching frequency F sw .

Fig. 9

 9 Fig. 9 describes a circuital model for actual inductors usually used for the circuit simulation tools. This model includes an ideal inductor L, a series resistance R s , a parasite capacitance C p , and a parallel losses resistance R p .

Figure 9 :

 9 Figure 9: Boost converter with inductor model.

Figure 10 :

 10 Figure 10: Impact of the SRF on the drain current. Color nomenclature: blue -drain current for SRF =250MHz, red -drain current for SRF =30MHz.

Fig. 11

 11 Fig. 11 depicts the frequency behavior of the associated inductor losses R f and the quality factor Q f for a power inductor of 8.2µH, Q(30MHz) = 100 and SRF = 250MHz. Fig. 12 shows the inductor current for the aforementioned specifications in comparison with an ideal inductor. The frequency distribution of power losses are plotted in Fig. 13. This figure shows that

Figure 11 :

 11 Figure 11: Frequency behavior of R f and Q f for the modeled inductor.

Figure 12 :

 12 Figure 12: Inductor current. Color nomenclature: blue -ideal inductor with low series resistance, red -power inductor considering Q(30MHz) = 100 and SRF = 250MHz.

Figure 13 :

 13 Figure 13: Inductor power distribution in the frequency-domain for Q(30MHz)=100 and SRF = 250MHz.

Figure 14 : 5 . 2 .

 1452 Figure 14: Real part of inductor impedance (R f ). Inductor reference SRP5015TA-8R2M.

Fig. 15

 15 Fig. 15 depicts the waveforms for square signals of case (a) for 10MHz

  Figure 15: Test inductor SRP5015A at 10MHz and 30MHz.

Figure 16 :

 16 Figure 16: FFT of inductor current for test at 10MHz -square signal. Inductor reference SRP5015A.

Figure 17 :

 17 Figure 17: Experimental results for 1MHz boost converter. Inductor SRR1210-8R2. P avg =492mW.

Figure 18 :

 18 Figure 18: FFT for the inductor current of 1MHz boost converter.

( a )Figure 20 :

 a20 Figure 19: Inductor current and voltage for boost converter 60V -10W at 30MHz

Figure 21 :

 21 Figure 21: Output voltage of boost converter 400V -400W at 30MHz.

Table 1 :

 1 Characterization and modeling of inductors

			Measure			Calculation
	Inductor Ref.	L	Rs	Q	SRF	Rp	Cp MAPE
		(µH) (mΩ) (10MHz)	(MHz) (KΩ) (pF)
	SRP5015TA-8R2M 7.8	190	25.5	24.2	15.3	5.6	0.28
	SRR1210-8R2Y	7.3	4.1	15.4	21.1	9.1	7.9	0.27
	7447713082	7.7	29	13.2	29.8	7.2	3.7	0.28
	744314850	8.2	24	7.7	35.4	4.3	2.5	0.15

Table 2 :

 2 Inductor power frequency-domain (SRP5015A). Test at 10MHz -square signal.

	Freq.(MHz)	0	10 30 50 70 90 Total
	P h (µW)	0.03 197 14 4.8 2.9 1.7 220

Table 3 :

 3 Inductor power in the time and frequency

		SRP5015TA-8R2M	SRR1210-8R2Y
		10MHz	30MHz	10MHz 30MHz
	Power time (µW)	239	848	257	925
	Power freq.(µW)	220	760	235	857
		7447713082	744314850
		10MHz	30MHz	10MHz 30MHz
	Power time (µW)	285	1380	299	818
	Power freq.(µW)	262	1140	283	767

Table 4 :

 4 Inductor power in the frequency domain

	Freq.(MHz) 0	1	2	3	4	5 Total
	P h (mW)	14 418 0.4 11 0.05 2.7 447
	5.4. Boost converter for 400V -400W at 30MHz

Acknowledgments

This work has been partially funded by the Region Occitanie Pyrenees-Mediterranee.


Joel Guiot

Michelle Leydet

• Guiot

J ( Cerege

Climatic variations across the Mediterranean

published or not. The documents may come


Maud Fagny

Frédéric Austerlitz

Understanding the adaptation of polygenic traits: the importance of gene regulatory networks

Keywords: Quantitative traits, Polygenic selection, Gene regulatory networks, Systems Biology, Adaptive evolution, Population genetics modelling Highlights

The adaptation of populations to local environments often relies on the selection of beneficial versions of polygenic traits. At the molecular level, this process corresponds to polygenic selection models, in which many independent loci, often located in regulatory regions, experience subtle shifts in frequency. Here, we first review the different approaches developed to detect signatures of polygenic selection, their power and limitations. We then highlight the contribution of systems biology to the understanding of the molecular bases of polygenic traits and the evolution of gene regulatory networks involved in these traits. Finally, we discuss the need for a unifying framework merging the fields of population genetics, quantitative genetics and systems biology to better understand the molecular bases of polygenic traits adaptation.

• The topology of gene regulatory elements put constraints on which genes and regulatory regions can be leveraged by polygenic selection.

Local adaptation of polygenic traits: the importance of regulatory regions

Many adaptive phenotypes -i.e. phenotypes involved in adaptation to local environments -are polygenic (see Box 1). Also called quantitative or complex traits, these phenotypes are determined by multiple genes and regulatory loci. This complex genetic architecture has strong implications on how these traits evolve. If these traits are under selection (either stabilising or directional), the underlying genetic model will be a model of polygenic selection, i.e. subtle shifts in frequencies at many independent loci coding for each adaptive trait [START_REF] Pritchard | Adaptation -not by sweeps alone[END_REF][START_REF] Fraser | Gene expression drives local adaptation in humans[END_REF][START_REF] Höllinger | Polygenic adaptation: From sweeps to subtle frequency shifts[END_REF]. This is for example the case for height, immune response, and metabolism in humans [START_REF] Hancock | Adaptations to climate-mediated selective pressures in humans[END_REF][START_REF] Daub | Evidence for polygenic adaptation to pathogens in the human genome[END_REF], flowering time, and stress response in plants, [START_REF] He | The Footprint of Polygenic Adaptation on Stress-Responsive Cis -Regulatory Divergence in the Arabidopsis Genus[END_REF][START_REF] Zan | A Polygenic Genetic Architecture of Flowering Time in the Worldwide Arabidopsis thaliana Population[END_REF][START_REF] Hämälä | Gene Expression Modularity Reveals Footprints of Polygenic Adaptation in Theobroma cacao[END_REF], or pathogenicity in yeast [START_REF] Fraser | Systematic Detection of Polygenic cis-Regulatory Evolution[END_REF]. Moreover, in humans, a lot of disorders such as autoimmune diseases and metabolic disorders are considered to be consequences of past adaptation of polygenic traits to changes in local environments, such as pathogens or lifestyle transition [START_REF] Barreiro | From evolutionary genetics to human immunology: how selection shapes host defence genes[END_REF][START_REF] Chouchani | Metabolic adaptation and maladaptation in adipose tissue[END_REF]. Being able to detect polygenic adaptation is thus important not only to understand the evolution of polygenic traits, but also the origin of many human diseases.

Our apprehending of the effects of directional selection on polygenic traits at the molecular level has strongly improved in the past ten years. QTL and genome-wide association studies (GWAS), have led to the identification of genetic variants associated with these polygenic traits, which were located overwhelmingly in non-coding regions [START_REF] Edwards | Beyond GWASs: illuminating the dark road from association to function[END_REF][START_REF] Rodgers-Melnick | Open chromatin reveals the functional maize genome[END_REF] -up to 90% in some species -highlighting the importance of regulatory regions. Population genetics studies, leveraging the ever growing number of genomic data sets available and the improvements in genomic functional annotation have confirmed that local adaptation often relies on selection of mutations located in regulatory regions that impact gene expression levels [START_REF] Fraser | Systematic Detection of Polygenic cis-Regulatory Evolution[END_REF][START_REF] Enard | Genome-wide signals of positive selection in human evolution[END_REF]. This has spurred efforts to develop population genetic approaches to detect signatures of selection beyond coding genomic regions. Independently, these results helped creating new theoretical models, in which polygenic traits result from a complex interplay of regulatory relationships between genes and different molecular actors at the cellular level. This led to the development of systems biology models aiming at understanding how these regulatory relationship evolve under directional selection, which brought new clues to understand the molecular bases of polygenic traits adaptation to local environments.

Detecting genomic signatures of polygenic selection: an ongoing quest

During the past few years, population geneticists have attempted to develop approaches to detect signatures of polygenic selection on genomic data. The first difficulty is to identify the type of signatures that one must look for. For this purpose, models have been developed to investigate the expected impact of polygenic adaptation on population genomic data [START_REF] Höllinger | Polygenic adaptation: From sweeps to subtle frequency shifts[END_REF][START_REF] Latta | Differentiation of allelic frequencies at quantitative trait loci affecting locally adaptive traits[END_REF][START_REF] Le Corre | Genetic variability at neutral markers, quantitative trait loci and trait in a subdivided population under selection[END_REF][START_REF] Kremer | Decoupling of differentiation between traits and their underlying genes in response to divergent selection[END_REF][START_REF] Chevin | Selective sweep at a quantitative trait locus in the presence of background genetic variation[END_REF][START_REF] Stephan | Signatures of positive selection: from selective sweeps at individual loci to subtle allele frequency changes in polygenic adaptation[END_REF] (see Box 2).

They assume additive and/or epistatic interactions among loci. These studies have shown mainly that polygenic adaptation is likely to yield subtle frequency shifts at the genes involved in the polygenic trait: polygenic adaptation events target several mutations simultaneously and leave weak molecular signatures on the genome [START_REF] Pritchard | The genetics of human adaptation: hard sweeps, soft sweeps, and polygenic adaptation[END_REF]. Recently, Hollinger and colleagues have however shown that the type of genomic signature observed ranges from sweeps (either total or partial) to subtle frequency shifts, depending on the population-scaled mutation rate [START_REF] Höllinger | Polygenic adaptation: From sweeps to subtle frequency shifts[END_REF].

Methods to detect genomic signatures of polygenic selection events have mainly focused on combining signatures of positive selection across sets of mutations that are defined a priori as likely candidates for adaptive pressures. In particular, some studies have proposed to group mutations for which allele frequencies correlate with specific ecological variables [START_REF] Hancock | Adaptations to climate-mediated selective pressures in humans[END_REF][START_REF] Haasl | Fifteen years of genomewide scans for selection: trends, lessons and unaddressed genetic sources of complication[END_REF]. These approaches, while powerful, require access to past and present environmental data for the populations under study, which are often difficult to obtain. Other approaches use gene pathways and information about biological functions or GWAS results to group mutations together [START_REF] Daub | Evidence for polygenic adaptation to pathogens in the human genome[END_REF][START_REF] Hämälä | Gene Expression Modularity Reveals Footprints of Polygenic Adaptation in Theobroma cacao[END_REF][START_REF] Berg | A Population Genetic Signal of Polygenic Adaptation[END_REF]. The first approach consists in investigating whether known pathways or gene sets are particularly enriched in selection signals [START_REF] Daub | Evidence for polygenic adaptation to pathogens in the human genome[END_REF], with many loci exhibiting a high level of genetic differentiation among populations. A recently developed related approach looks for local subnetworks under strong differentiation [START_REF] Gouy | Detecting gene subnetworks under selection in biological pathways[END_REF].

The second approach consists in testing whether a known trait has been submitted to divergent selection [START_REF] Berg | A Population Genetic Signal of Polygenic Adaptation[END_REF][START_REF] Turchin | Evidence of widespread selection on standing variation in Europe at height-associated SNPs[END_REF][START_REF] Racimo | Detecting polygenic adaptation in admixture graphs[END_REF][START_REF] Robinson | Population genetic differentiation of height and body mass index across Europe[END_REF]. These methods require that previously performed genome-wide association studies (GWAS) are available, which have detected loci associated with the traits and estimated their additive effect on these traits. This allows computing the polygenic score of the genotyped individuals in the studied populations, which can then be compared to what is expected under genetic drift only. Because they rely on already existing GWAS studies, these methods cannot be used blindly on the genome.

On the other hand, other methods have focused on scanning the genome in order to detect loci with high level of differences in allele frequencies among populations, without a priori grouping mutations. This is the case for example of the SDS approach [START_REF] Field | Detection of human adaptation during the past 2000 years[END_REF], OutFLANK [START_REF] Whitlock | Reliable detection of loci responsible for local adaptation: Inference of a null model through trimming the distribution of F st[END_REF] and PCAdapt [START_REF] Duforet-Frebourg | Detecting Genomic Signatures of Natural Selection with Principal Component Analysis: Application to the 1000 Genomes Data[END_REF]. These genome-scan methods can be used blindly on the genome, however they were mainly developed in a context of strong selection on a single locus and simulation studies have shown that their efficiency for detecting loci evolving under polygenic selection is rather limited [START_REF] De Villemereuil | Genome scan methods against more complex models: when and how much should we trust them?[END_REF][START_REF] Vilas | A simulation study on the performance of differentiation-based methods to detect selected loci using linked neutral markers[END_REF].

Shortcomings of polygenic selection detection approaches

The shortcomings of existing approaches to detect polygenic selection have been extensively discussed in recent papers [START_REF] Novembre | Tread Lightly Interpreting Polygenic Tests of Selection[END_REF][START_REF] Rosenberg | Interpreting polygenic scores, polygenic adaptation, and human phenotypic differences[END_REF][START_REF] Berg | Reduced signal for polygenic adaptation of height in UK Biobank[END_REF][START_REF] Sohail | Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies[END_REF] that have raised several caveats concerning the results of these polygenic selection tests. The first set of caveats are technical and relate primarily to GWAS-based results: if not correctly taken into account, a hidden population structure in the samples used for GWAS can lead to the false detection of directional polygenic selection signal on the trait under study. This has been demonstrated for height in European human populations [START_REF] Berg | Reduced signal for polygenic adaptation of height in UK Biobank[END_REF][START_REF] Sohail | Polygenic adaptation on height is overestimated due to uncorrected stratification in genome-wide association studies[END_REF]. Similarly, GWAS ascertainment biases can confuse selection signals [START_REF] Uricchio | Evolutionary perspectives on polygenic selection, missing heritability, and GWAS[END_REF]. Transferring GWAS results across time and space can also lead to misidentification of causal mutations, the impact of which on the selection tests is uncertain.

The last set of caveats relates to the interpretation of test results. They include the potential existence of gene x environment interactions and the poorly known interactions between traits [START_REF] Novembre | Tread Lightly Interpreting Polygenic Tests of Selection[END_REF], both of which can hide the selection signals or inverse the direction of selection. Finally, an important part of the selection signals might be lost in current studies because a fundamental component of the genetic architecture of polygenic traits is often minimized or totally ignored: epistatic interactions [START_REF] Csillery | Detecting the genomic signal of polygenic adaptation and the role of epistasis in evolution[END_REF]. On this topic, Payne and colleagues recently highlighted the importance of studying the topography of the adaptive landscape and, in particular, the exceptional characteristics of gene expression regulatory networks in terms of epistasic interactions, pleiotropy and robustness [START_REF] Payne | The causes of evolvability and their evolution[END_REF].

To this list of caveats, we can add another one: while carrying the majority of the genetic variation associated with polygenic traits, the regulatory regions are a blind spot of the methods designed to detect polygenic selection described above. These approaches, which combine signatures of positive selection across loci associated to a given trait, often only target genes, completely ignoring regulatory regions. The recently developed approaches based on GWAS results include some genetic variants from regulatory regions. However, because they only include the variants that are statistically significantly associated to the phenotype at the genome-wide level, they miss a good part of the genetic bases of polygenic phenotypes -the famous "missing heritability" [START_REF] Manolio | Finding the missing heritability of complex diseases[END_REF][START_REF] Boyle | An Expanded View of Complex Traits: From Polygenic to Omnigenic[END_REF].

On the other hand, blind genome-wide scans have detected the signatures of polygenic selection in regulatory regions. Several studies have indeed revealed an enrichment in positive selection signals among cis-regulatory elements and expression QTLs [START_REF] Kudaravalli | Gene expression levels are a target of recent natural selection in the human genome[END_REF][START_REF] Torgerson | Evolutionary processes acting on candidate cis-regulatory regions in humans inferred from patterns of polymorphism and divergence[END_REF][START_REF] Rougeux | Polygenic selection drives the evolution of convergent transcriptomic landscapes across continents within a Nearctic sister species complex[END_REF]. However, the lack of functional annotation of these regions often prevents identification of targeted phenotypes, beyond general results of Gene Ontology enrichment assays.

In short, despite recent advances in polygenic selection signature testing, several difficulties persist. While the technical caveats point to the need to carefully consider potential confounding effects and to assess the sensitivity and specificity of polygenic selection tests, most of the difficulties in interpreting polygenic selection tests can be explained by the "genomic only" approaches used in these tests. These difficulties indeed arise from our limited knowledge of the true genetic architecture and molecular bases of polygenic traits. These traits are determined by various regulatory regions, and depend on a complex interplay between different actors that are expressed following a rigid spatio-temporal framework. Ignoring these features makes it difficult to interpret results of polygenic selection tests, and to identify the biological functions and traits potentially under polygenic selection.

Embracing the Complexity: Contributions from Systems Biology

Organizing and integrating information from diverse sources and understanding how polygenic phenotypes are generated from the genome at different levels, from the cell to the organism, is the core of systems biology [START_REF] Kirschner | The Meaning of Systems Biology[END_REF]. The rise of this research field initiated a switch in the perception of the molecular bases of polygenic traits, from a gene-first to an interaction-first model. Quantitative genomics apprehends indeed the molecular bases of polygenic phenotypes as a collection of independent genes, each coding for a fraction of the phenotype, with potentially some interactions.

Conversely, systems biology focuses primarily on these interactions, modelling the regulatory relationship between genes, their regulators and their products. During the past twenty years, buoyed by the availability of high-throughput data, numerous approaches have been developed to model these regulatory relationships and to efficiently deal with their complexity. In particular, the representation of these relationships as Gene Regulatory Networks (GRNs), which gather tens to thousands of genes, has allowed better understanding the molecular bases of polygenic traits [START_REF] Sonawane | Network Medicine in the Age of Biomedical Big Data[END_REF].

From the systems biology point of view, understanding the molecular bases of polygenic selection consists, therefore, in deciphering how regulatory interactions among genes evolve under constraints. Models have been developed to study how GRNs evolve, which explicitly simulates these interactions. Initially developed by Wagner, they model gene networks through a matrix whose coefficients determines how each gene regulates each other and itself [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF][START_REF] Wagner | Does evolutionary plasticity evolve?[END_REF]. These coefficients vary among individuals and can evolve through mutations that change the regulatory relationships between genes. Each matrix leads to a given level of expression of each locus. Selection will occur through the assumption that there is an optimal level of expression for each gene.

This model was first developed to study the evolution of gene duplications [START_REF] Wagner | Evolution of gene networks by gene duplications: a mathematical model and its implications on genome organization[END_REF].

It allowed then discovering some properties of GRN evolution. One of the major results is the demonstration that, under stabilizing selection, network evolves toward a state where a single mutation becomes less likely to reduce the fitness of the individuals [START_REF] Wagner | Does evolutionary plasticity evolve?[END_REF], a phenomenon called canalisation [START_REF] Waddington | Canalization of development and the inheritance of acquired characters[END_REF]. This insensitivity to mutations becomes stronger when the number of connections in the network increases [START_REF] Siegal | Waddington's canalization revisited: Developmental stability and evolution[END_REF] or when phenotypic optima become more extreme [START_REF] Runneburger | Why and how genetic canalization evolves in gene regulatory networks[END_REF]. Conversely, this robustness of the network promotes the accumulation of hidden loss-of-function mutations, which can then accelerate adaptation in a changing environment [START_REF] Bergman | Evolutionary capacitance as a general feature of complex gene networks[END_REF].

GRN evolution models have also helped better understanding how the structure of these GRNs impacts the individual response of each gene to selection, as a function of selection intensity and distance between the current phenotype and the optimal one. Genes that are more upstream in the network (i.e. strongly regulating the other genes and/or less regulated by them) are more affected by selection [START_REF] Rhone | Impact of selection on genes involved in regulatory network: a modelling study[END_REF]. This could explain the empirical observation that transcription factors are evolving faster than the other genes [START_REF] Bustamante | Natural selection on protein-coding genes in the human genome[END_REF].

Impact of GRN structure on architecture of polygenic selection

Understanding the evolution of GRNs that underlies polygenic phenotypes is currently a hot topic in evolutionary biology. The various scenarii of GRN evolution observed have been summarized for animals, insects [START_REF] Halfon | Perspectives on Gene Regulatory Network Evolution[END_REF] and plants [START_REF] Das Gupta | Gene networks and the evolution of plant morphology[END_REF]. Modifications of GRNs seems to be frequent during speciation [START_REF] Ten Tusscher | The role of genome and gene regulatory network canalization in the evolution of multi-trait polymorphisms and sympatric speciation[END_REF] and domestication [START_REF] Doebley | The Molecular Genetics of Crop Domestication[END_REF][START_REF] Swanson-Wagner | Reshaping of the maize transcriptome by domestication[END_REF], and can lead to dramatic changes in developmental patterns and other polygenic phenotypes. Moreover, there are many examples of directional and polygenic selection events affecting the gene regulatory networks [START_REF] Hämälä | Gene Expression Modularity Reveals Footprints of Polygenic Adaptation in Theobroma cacao[END_REF][START_REF] Fraser | Systematic Detection of Polygenic cis-Regulatory Evolution[END_REF][START_REF] Rougeux | Polygenic selection drives the evolution of convergent transcriptomic landscapes across continents within a Nearctic sister species complex[END_REF][START_REF] Bullard | Polygenic and directional regulatory evolution across pathways in Saccharomyces[END_REF]]. The analysis of cellular GRN and their evolution has led to hypotheses about how polygenic selection can act, and how gene regulatory network structure might facilitate or constrain it. Three GRN features are worthy of attention in this context: redundancy, modularity, and their scale-free property (see Figure 1 and Box 3).

Redundancy is a characteristics of GRNs that ensures stability of gene expression, buffers noise and ensures phenotype robustness to random mutations that may cause gene inactivation [START_REF] Runneburger | Why and how genetic canalization evolves in gene regulatory networks[END_REF][START_REF] Macneil | Gene regulatory networks and the role of robustness and stochasticity in the control of gene expression[END_REF].

It has been shown also to facilitate polygenic adaptation. For example, an evolution experiment in Drosophila showed phenotypic convergence in response to a temperature change in ten populations, while the genetic analysis revealed that different pathways were targeted by polygenic selection in each replicate [START_REF] Barghi | Genetic redundancy fuels polygenic adaptation in Drosophila[END_REF]. Redundancy of GRNs thus facilitated polygenic selection by increasing the number of available targets to reach the optimum phenotype (see Figure 1B). However, this polygenic selection will be difficult to detect, as it will potentially affect different pathways among populations confronted to the same constraints.

GRNs are also highly modular, i.e., made of small groups of genes strongly interconnected with each other but only loosely connected to the rest of the network, thus constituting a module. It has been proposed, in the omnigenic model, that as all genes are seemingly interconnected within the GRN, each trait might be influenced by all genes [START_REF] Boyle | An Expanded View of Complex Traits: From Polygenic to Omnigenic[END_REF]. Theoretically, such extreme generalized pleiotropy -property of a gene that influences several, apparently unrelated traits -would severely constrain the possibility for polygenic selection to act, because it could not modify a given trait without influencing all the others. However, the modularity of GRNs can mitigate this effect. Indeed, in each module, only a few genes are connected with genes from other modules. Most genes in a module are thus isolated from the rest of the network and can be leveraged by polygenic selection (see Figure 1A). Finally, GRNs tend to be scale-free, i.e. there are many more genes that are connected to only a few other genes than there are genes connected to lots of other genes [START_REF] Albert | Scale-free networks in cell biology[END_REF] (see Figure 1B). This has been shown on various types of networks [START_REF] Sonawane | Network Medicine in the Age of Biomedical Big Data[END_REF][START_REF] Albert | Scale-free networks in cell biology[END_REF][START_REF] Platig | Bipartite Community Structure of eQTLs[END_REF]. As confirmed by observations on molecular data (transcriptomics or proteomics), the level of connection of a gene in a network is inversely correlated with its level of genetic diversity [START_REF] Jordan | Conservation and Coevolution in the Scale-Free Human Gene Coexpression Network[END_REF][START_REF] Josephs | The Relationship between Selection, Network Connectivity, and Regulatory Variation within a Population of Capsella grandiflora[END_REF][START_REF] Luisi | Recent Positive Selection Has Acted On Genes Encoding Proteins With More Interactions Within the Whole Human Interactome[END_REF]. This means that the handful of highly pleiotropic genes at the center of these networks often evolve under strong negative selection and present few variations of expression among individuals. Conversely, the vast majority of genes are at the periphery of the network, involved in only one or few pathways. Changes in their regulation will, therefore, only very slightly affect the other pathways (see Figure 1A andB). This leaves open a large field full of candidates for polygenic selection to act on, both at the level of the gene sequence and the gene regulatory elements, while limiting the potential side effects of selection on other pathways.

Polygenic selection and GRNs: advantages and limits

Theoretical models of GRNs combined with results on real data sets allow better understanding the molecular bases of polygenic phenotypes and their evolution. They also help understanding where and how polygenic selection can act, and where to look for genomic signatures of selection.

However, the models often focus on the evolution of unipartite networks of genes interacting with each other, thus overlooking important features that would be useful to get a complete picture of how polygenic selection can impact GRNs. First, despite the evidences for pervasive polygenic selection targeting regulatory mutations [START_REF] Fraser | Systematic Detection of Polygenic cis-Regulatory Evolution[END_REF], theses models do not include explicitly how modifications in non-coding regulatory regions, such as enhancers, can impact the expression level of several genes at once, while not modifying their interactions. An important source of polygenic trait variation and evolution is thus not taken into account. Some current approaches integrate heterogeneous data using multipartite networks, in order to identify not only gene-gene regulatory interactions but also interactions between regulatory regions and genes [START_REF] Sonawane | Network Medicine in the Age of Biomedical Big Data[END_REF][START_REF] Fagny | Exploring regulation in tissues with eQTL networks[END_REF][START_REF] Baker | Integrating transcriptomic network reconstruction and eQTL analyses reveals mechanistic connections between genomic architecture and Brassica rapa development[END_REF].

These approaches could serve as bases for the development of new GRN evolution models. Second, most studies on GRN evolution do not model individual nucleotides. Consequently, they allow to study questions such as redundancy, pleiotropy and canalisation, but not to dissect the genetic bases of evolution, and to understand the genomic signatures that may be expected in different polygenic selection scenarios.

Concluding Remarks and Future Perspectives: towards systems population genetics

In order to detect polygenic selection and to disentangle the impact of several confounding factors, we need to better understand not only the genetic architecture of polygenic traits, but also their molecular architecture. A new integrative evolutionary framework becomes quite indispensable, which would take this architecture into account, in particular its redundancy and the fact that two selective events in the same population might not lead to the same selection scenario at the molecular level [START_REF] Barghi | Polygenic adaptation: a unifying framework to understand positive selection[END_REF]. Here, we argue that, in order to achieve such a unifying framework, we need not only to identify which genes and pathways contribute to polygenic traits, and how redundant they are, but also investigate the complex regulatory interactions that govern their expression (see Outstanding Questions). In particular, we need to take into account how the genomic regulatory regions impact the expression of these traits. We should focus on the crucial post-trancriptional regulatory steps, on the gene-gene interactions, and maybe more importantly, on the position of these genes in the cellular gene regulatory network. All of this could be achieved by combining population genomics and system biology approaches.

Outstanding Questions

• Redundancy allows for the selection of the same version of a trait through the selection of genetic variants in different pathways. Do we observe this phenomenon at the scale of the population? In that case, which molecular signatures do such evolutionary processes leaves on the genome and how can we detect them?

• Pleiotropy is a major property of gene regulatory networks actors that constrains how polygenic selection acts. Genes are often involved in several distinct pathways. What is the level of pleiotropy of the regulatory elements, in particular of the distal ones such as enhancers and how does it constrain their evolution?

• Classical population genetics approaches model the evolution of genetic sequences but ignore their complex interactions in gene regulatory networks. On the other hand, gene regulatory networks model the evolution of these regulatory relationships but not their genetic bases.

Will more integrated approaches that explicitly model the molecular evolution of regulatory regions in the context of gene regulatory networks improve our understanding of the molecular bases of local adaptation?

• Selection on polygenic traits yields linkage disequilibrium among genes involved in the traits, at the within-and among-population levels [START_REF] Latta | Differentiation of allelic frequencies at quantitative trait loci affecting locally adaptive traits[END_REF]. How do regulatory interactions affect this linkage disequilibrium? Can it be used to detect polygenic selection?

• GWAS-based method aiming at detecting polygenic selection use the additive effect of each gene on the trait. Could the interaction (dominance, epistasy) between these genes also be integrated in these approaches? Box 1. Genetic architecture of traits: from the Mendelian to the omnigenic model Polygenic traits have been the focus of many theoretical and empirical studies since Fisher's seminal model of polygenic inheritance, which bridged the gap between the Mendelian inheritance model and the observation of continuous traits. This model suggests that the apparent continuity of many phenotypic traits can be explained if multiple genes contribute additively to these traits [START_REF] Fisher | XV.-The Correlation between Relatives on the Supposition of Mendelian Inheritance[END_REF]. Classic genetic studies, recent advances in quantitative genetics and genome-wide association studies (GWAS) have demonstrated the existence of a continuum of genetic architectures for traits, in many organisms including yeast, insects, worms, plants and mammals [START_REF] Manolio | Finding the missing heritability of complex diseases[END_REF][START_REF] Mackay | The genetics of quantitative traits: challenges and prospects[END_REF][START_REF] Holland | Genetic architecture of complex traits in plants[END_REF][START_REF] Timpson | Genetic architecture: the shape of the genetic contribution to human traits and disease[END_REF]. These architectures range from the Mendelian model, where a trait is determined by a single gene (e.g. ABO blood groups in primates, Mendel's peas color and shape) to the infinitesimal model [START_REF] Barton | The infinitesimal model: Definition, derivation, and implications[END_REF] where a seemingly infinite number of loci are involved in the determinism of a trait (e.g. height, in both humans and plants [START_REF] Turchin | Evidence of widespread selection on standing variation in Europe at height-associated SNPs[END_REF][START_REF] Peiffer | The Genetic Architecture Of Maize Height[END_REF]). This last concept has been highlighted in a recent paper [START_REF] Boyle | An Expanded View of Complex Traits: From Polygenic to Omnigenic[END_REF] that labelled it the 'omnigenic' model. In between, we find a myriad of traits that present an oligogenic or polygenic determinism, such as flowering time in plants, tolerance to ethanol in yeasts or susceptibility to diabetes in humans [START_REF] Peiffer | The Genetic Architecture Of Maize Height[END_REF][START_REF] Pais | Comparative Polygenic Analysis of Maximal Ethanol Accumulation Capacity and Tolerance to High Ethanol Levels of Cell Proliferation in Yeast[END_REF][START_REF] Salome | Genetic Architecture of Flowering-Time Variation in Arabidopsis thaliana[END_REF][START_REF] Fuchsberger | The genetic architecture of type 2 diabetes[END_REF].

Box 2. Polygenic selection

Polygenic selection consists in natural or artificial selection occurring on a trait coded by several loci. Under this framework, even when the phenotypic optimum differs strongly from the current average phenotype, populations can evolve toward this optimum without substantial changes in allelic frequencies [START_REF] Latta | Differentiation of allelic frequencies at quantitative trait loci affecting locally adaptive traits[END_REF]. Indeed, each locus can undergo only swift changes in these frequencies, while allelic associations will appear among the different loci, allowing to reach the phenotypic optimum. In this case, methods aiming at detecting signals of selective sweeps such as iHS [START_REF] Voight | A map of recent positive selection in the human genome[END_REF] or nSL [START_REF] Ferrer-Admetlla | On detecting incomplete soft or hard selective sweeps using haplotype structure[END_REF] will probably not detect any signal, as no locus is undergoing a strong shift in allelic frequencies. Note, however, that this will depend on the number of loci involved in the trait: if this number is low, each locus will respond more strongly in terms of allelic frequencies changes and some might be detected. Moreover, if among the loci coding for the trait, some are major loci with a high contribution to genetic variance, these loci are more likely to be affected by a strong shift in allelic frequencies [START_REF] Le Corre | Genetic variability at neutral markers, quantitative trait loci and trait in a subdivided population under selection[END_REF][START_REF] Kremer | Decoupling of differentiation between traits and their underlying genes in response to divergent selection[END_REF], and thus to be detected by selective sweeps tests.

Other approaches modeled the fate of a locus coding for a given trait, while all other loci coding for this trait were treated as background genetic variance [START_REF] Chevin | Selective sweep at a quantitative trait locus in the presence of background genetic variation[END_REF]. They find that a selective sweep will be observed in some cases at the focal locus, depending upon its relative variance as compared to the background loci and also upon the shape of the selection function. In a model that considered also epistatic interactions, Hollinger and colleagues then showed that the probability of observing a sweep at a locus involved in a quantitative trait will depend mainly on a parameter called the population-scaled background mutation rate Θ bg , i.e. the product of the effective size of the population by the total mutation rate of all background loci [START_REF] Höllinger | Polygenic adaptation: From sweeps to subtle frequency shifts[END_REF]. Low values for Θ bg will lead to the occurrence of selective sweeps at the focal locus, while high values will yield small polygenic allele frequency shifts. Several partial selective sweeps will be observed for intermediate values.

Box 3. Gene Regulatory Networks

Gene Regulatory Networks (GRNs) describe a collection of molecules that interact together to regulate the level of expression of genes, both as mRNA and proteins. These regulatory relationships can be direct or indirect, for example when a transcription factor increases another protein level through the activation of its coding gene transcription.

GRNs are often represented as graphs, mathematical objects that represent the elements (nodes) and the interactions between elements (edges). They can include one type of elements such as proteins (unipartite graphs), or several types of elements such as transcription factors, regulatory mutations and mRNAs (multipartite graphs). The study of the topology of these graphs allows better understanding how the elements interact together and how GRNs evolve. GRNs, whether in yeasts, insects, plants or mammals, show some shared topological characteristics. Most importantly, they are scale-free. This means that very few of the nodes within the network are connected to many others -the key or master regulators -while most of the nodes are connected only to a few other nodes (see Figure Figure 1B). GRNs are also highly modular, i.e. organized into groups of nodes that are strongly interconnected but only loosely connected to the rest of the network, thus constituting a module. The modules often correspond to biological pathways.

The topology of the network also seems to constrain the evolutionary history of each node. In GRNs, the handful of high-degree genes, i.e. connected to many other genes, are often known to be pleiotropic, involved in many pathways, and to be very constrained [START_REF] Jordan | Conservation and Coevolution in the Scale-Free Human Gene Coexpression Network[END_REF], evolve under strong negative selection and present few variations of expression among individuals. Conversely, peripherical genes, which are only connected to few other genes, are often included in only one pathway.

They are also more likely to evolve under positive selection. products. Genes are colored by topological characteristics: high degree genes in red, high core-score genes in yellow and others in white. The network is structured into modules, i.e. groups of genes that regulate each others. These modules often correspond to biological functions. Topological characteristics of gene regulatory networks such as redundancy (several modules involved in the regulation of the same biological functions) and pleiotropy (genes involved in the regulation of several biological functions) are highlighted. (B) Expected power law distribution of the node degree in q gene regulatory network.

Gene Expression

Regulatory regions

Bipartite network

Integrated approach: Identify regulatory modules evolving under polygenic selection

Selected regulatory module

Population and quantitative genomics:


Louis Duchemin

Mathilde Paris

Bastien Boussau
email: bastien.boussau@univ-lyon1.fr

Predicted effects of summer holidays and seasonality on the SARS-Cov-2 epidemic in France

The SARS-CoV-2 epidemic in France has had a large death toll. It has not affected all regions similarly, since the death rate can vary several folds between regions where the epidemic has remained at a low level and regions where it got an early burst. The epidemic has been slowed down by a lockdown that lasted for almost eight weeks, and individuals can now move between metropolitan French regions without restriction. In this report we investigate the effect on the epidemic of summer holidays, during which millions of individuals will move between French regions. Additionally, we evaluate the effect of strong or weak seasonality and of several values for the reproduction number on the epidemic, in particular on the timing, the height and the spread of a second wave. To do so, we extend a SEIR model to simulate the effect of summer migrations between regions on the number and distribution of new infections. We find that the model predicts little effect of summer migrations on the epidemic, because the number of migrating infectious individuals are low as a consequence of the lockdown. However, all the reproduction numbers above 1.0 and the seasonality parameters we tried result in a second epidemic wave, with a peak date that can vary between October 2020 and April 2021. If the sanitary measures currently in place manage to keep the reproduction number below 1.0, the second wave will be avoided. If they keep the reproduction number at a low value, for instance at 1.1 as in one of our simulations, the second wave is flattened and could be similar to the first wave.

Introduction

On June 29 2020, the World Health Organization documented that more than 10 million cases of infection by SARS-Cov-2 had been reported (World Health Organization, 2020). In continental France, the epidemic seemed to be contained after a series of measures, including a drastic lockdown from March 17 to May 11. During the lockdown, one required a self-authorisation to leave home. The lockdown has successfully slowed down the spread of the epidemic, reducing the reproduction number R t from around 3 to less than 1 [START_REF] Salje | Estimating the burden of SARS-CoV-2 in France[END_REF][START_REF] Sofonea | Epidemiological monitoring and control perspectives: application of a parsimonious modelling framework to the COVID-19 dynamics in France. medRxiv[END_REF][START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF]. The lockdown was then lifted progressively. On May 11, workers were allowed to go back to work, and schoolchildren to go back to school, with new sanitary measures enforced in shops, at work and in schools. Individuals were not allowed to change 2 region or to move more than 100km away from their home. On June 22, all children were expected to attend school, and some of the sanitary measures were no longer enforced. In particular, restrictions on travel within continental France were lifted. As of June 25, the epidemic in continental France is considered under control by Santé Publique France (SPF), the governmental agency in charge of monitoring the epidemic (Santé Publique France, 2020). SPF estimates there are 35 new cases per day per 1 000 000 inhabitants. However, the epidemic has not had the same impact on all regions, and as a result the epidemic still does not have the same intensity in all regions. This can be seen in the predicted proportion of infectious individuals per region on July 1st 2020 (Fig. 1). These predictions are based on a Bayesian model that uses the number of SARS-CoV-2 related deaths through time to monitor the epidemic [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF]. This model was found to have good accuracy and agreed with other models of the epidemic in France. SPF also estimates a reproductive number at 0.92 (95% confidence interval: 0.89-0.96) for continental France for the week preceding June 25. Such a low R 0 is due to sanitary measures that are still in place, such as the prohibition of large gatherings or advice on hand washing and mask wearing in shops and at work. But it may also be due 3 to SARS-CoV-2 susceptibility to weather conditions [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF][START_REF] Kissler | Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period[END_REF], i.e. its seasonality. Other coronaviruses are known to be seasonal, with a higher reproduction number in winter than in summer, but it is still unknown how seasonal SARS-CoV-2 may be. In particular, in their investigations of SARS-CoV-2 seasonality, [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF] used a seasonality value of 0.5 for Northern temperate countries, while [START_REF] Kissler | Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period[END_REF] used a seasonality value of 0.13. Thus it is important to assess how these different estimates affect future trajectories for the epidemic in France.

Between 2014 and 2016, on average 119 million people moved between regions during holidays (Direction Générale des Entreprises, 2018) (Fig. 2). Most of these migrations happen during summer, and could help spread the virus between regions, including moving the virus into regions that comparatively had so far been spared by it. In this manuscript, we studied the effect summer holiday migrations between regions may have on the spread of the SARS-CoV-2 epidemic in France. We implemented a compartmental Susceptible-Exposed-Infectious-Recovered/Removed (SEIR) model to simulate the infectious process, and complemented it with migrations between regions as gathered from (Direction Générale des Entreprises, 2018). We initialized the number of individuals in each compartment based on simulations from [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF])'s 4 model. We first investigated how seasonality and reproduction numbers interact to determine the size and the timing of a second wave of infections in France. To this end, we examined two seasonality values and three reproductive number values for the summer. We finally investigated how migrations affected the number of infections in each region.

3 Material and methods

Models

SEIR model

We use a Susceptible-Exposed-Infectious-Recovered/Removed (SEIR) model to simulate the SARS-Cov-2 epidemic in France. Such models are used widely to simulate epidemics, and a formal presentation can be found in [START_REF] Hethcote | The Mathematics of Infectious Diseases[END_REF] or [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF]. A brief intuitive presentation follows. "Susceptible" individuals may become infected because they coexist with infectious individuals; upon infection, an individual becomes "Exposed". Exposed individuals then become "Infectious" after an incubation period. Infectious individuals contribute to the spread of the epidemic, and become "Recovered" or "Removed" after an infectious period. Individuals in the R category cannot move to another category.

The starting code base was taken from [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF], who used their SEIR model to investigate the influence of migrations and seasonality on viral spread worldwide. In their model, seasonality and migration parameters were drawn from prior distributions, and other parameters were fixed to values drawn from the literature. They used this model to examine possible trajectories that the worldwide epidemic could take over the next months and years. We built upon this model to investigate the influence of seasonality and summer holiday migrations on the epidemic in France. In particular, we used empirical data to inform the pattern of migrations between French regions in the mainland, and we take into account migrations in both directions: to the holiday region and back. We also used estimated reproduction numbers and combined them with a model of seasonality.

Simulating holiday migrations in the SEIR model

The SEIR model is based on differential equations that describe how individuals move between the S, E, I, R categories in a small amount of time dt. Each region has its own set of population size and starting numbers of individuals in the S, E, I, R categories. Other parameters such as the reproduction number or the infectious period are shared between regions. Contrary to agent-based models, each individual in the simulation is not considered separately from other individuals: instead, only the numbers of individuals in each compartment are considered. Simulation is performed by dividing the year in n intervals of time dt and running the differential equations repeatedly over 5 each interval. Migrations between regions happen every week during the summer holidays. A number of emigrant individuals are sent to destination regions and removed from their origin region; during the move, they keep their S, E, I, R category. These immigrants are added to the numbers of S, E, I, R individuals in the destination region and removed from the numbers of S, E, I, R individuals in the origin region for a week. Then, during a week, these immigrants can change category according to the process happening in the destination region. At the end of the week, the immigrants are brought back to their origin region. The total number of individuals coming back equals the number of individuals that have gone on holidays, but the proportions of S, E, I, R individuals among them can differ from the starting number. With such a process, individuals coming from regions with low infection rates can become infected if they move to regions with high infection rates, and infectious individuals migrating to another region may infect individuals in the destination region.

Joint simulation of seasonality and interventions

To simulate jointly non-pharmaceutical interventions and seasonality, we used estimates of R t from Duchemin et al.'s Bayesian model and modelled seasonality deterministically. Duchemin et al.'s Bayesian model uses the number of deaths in 13 regions in metropolitan France to estimate an R 0 value before the lockdown, and a set of other parameters that act as modifiers of the initial R 0 to provide estimates of daily R t during the lockdown, and after the lockdown. These R t values account for three elements: the effect of seasonality, the effect of other factors, notably non-pharmaceutical interventions, and the decline in the number of susceptible individuals. Taking the average of the R t values over a time period (i.e. before, during, or after the lockdown) provides an estimate that we will refer to as R average here and that factors out the effect of the decline in the number of susceptible individuals. But such an R average value still combines the effect of seasonality and the effect of non-pharmacetical interventions.

The following explains how we disentangle these two factors. We modelled seasonality deterministically as in [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF] as follows:

β(t) = β 0 × (1 + cos(2π(t -θ))) (1)
Here β(t) is the transmissibility at day t, i.e. the rate at which an infected individual infects another individual, is the strength of seasonal forcing, and θ is the date of peak transmissibility. We used the parameter β 0 , which was defined as the yearly average transmissibility in [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF] to separate the effects of seasonality and of other factors. To this end, we estimate one β average value per period, where each β average value corresponds to an R average value as follows:

R average = rec × β average (2)
Where rec is the recovery rate, i.e. the inverse of the time for an infected individual 6 to recover. Therefore, to estimate β average values, we use:

β average = rec × R average t i+1 t i (1 + cos(2π(t -θ)))dt (3)
Where t i and t i+1 are the boundary dates for period i. The resulting daily R t values are shown Fig. 3 for the period between July 1st 2020 and June 1st 2021. These were used as input for the SEIR model, which then handles the decline in the number of susceptible individuals on its own.

Data

Parameters

Since [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF]'s work, some parameters of the Covid19 disease have been determined more precisely, but others are still uncertain. In the case of reproduction numbers and seasonality, we compared several values to account for this uncertainty. We acknowledge a certain amount of arbitrariness in selecting the other parameter values but observed little sensitivity of the results when we introduced small variations.

We set the reproduction numbers by combining model-based estimates and seasonality estimates. We use the model-based estimates to provide average reproduction numbers R average over time-periods, and use seasonality estimates to allow variation 7 around these average values. We explored 3 different R average values for the period extending from July 1st 2020 to the end of the simulation, June 1st 2021: 1.1, 1.5 and 2.0. We did not explore values below 1.0 as they would be certain to lead to an extinction of the epidemic.

We considered two values for the seasonality of SARS-CoV-2. First we used a value of in agreement with [START_REF] Neher | Potential impact of seasonal forcing on a SARS-CoV-2 pandemic[END_REF]'s analyses of the seasonality of other coronaviruses in Northern temperate countries. Second we used a lower value of 0.13, in agreement with [START_REF] Kissler | Projecting the transmission dynamics of SARS-CoV-2 through the postpandemic period[END_REF]. In both cases, we considered that the maximum infectivity of the virus was on January 1st.

We set the infectious time, the time it takes between infection and becoming infectious, to 3 days, assuming it takes ≈ 5 days from infection to symptom onset [START_REF] Lauer | The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application[END_REF], and infectiousness starts ≈ 2 days before symptom onset [START_REF] He | Temporal dynamics in viral shedding and transmissibility of COVID-19[END_REF]. This is the time an individual remains in the compartment.

We set the recovery time, the time it takes between symptom onset and the end of infectiousness, to be 6 days, based on the inferred infection profile in [START_REF] He | Temporal dynamics in viral shedding and transmissibility of COVID-19[END_REF]. This is the time an individual remains in the compartment.

Initialization

Region wise proportions of S, E, I, R individuals on July 1 2020 were simulated using a Bayesian model of the French epidemic [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF]. We chose to use the model without mixture between regions as it provided predictions that were very similar to the mixture model and was easier to fit and simulate. This model was fitted with data up to June 29. The parameter values of each sample of the MCMC chain were used to simulate the epidemic spread up to July 1st 2020. Median values were extracted to initialize the SEIR model in each region. Exposed individuals were defined as individuals that had been infected less than 3 days ago. Infectious individuals were defined as individuals that had been infected less than 9 days ago but more than 3 days ago. Recovered/Removed individuals were defined as individuals that had been infected more than 9 days ago. All the other individuals were assigned to the Susceptible category. This assignment of individuals to categories is consistent with the SEIR parameters introduced above. The proportions of infectious individuals on July 1 2020 estimated by the model range between 3.1 × 10 -8 for Corsica and 6.9 × 10 -5 for Île de France, with a mean at 1.9 × 10 -05 . These numbers are in the same range as the number of new cases per day estimated by SPF in the week preceding June 25, which was 3.5 × 10 -5 .

Validation

Implementation and availability

The simulation is implemented in Python and is available at https://gitlab. in2p3.fr/boussau/futurecorona 4 Results

Differences between simulated numbers of infectious in-

dividuals using the SEIR model and estimates from a Bayesian model Fig. 4 shows that the SEIR model, when initialized with values drawn from [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF]'s model on March 1 2020, qualitatively captures the differences between regions that the Bayesian model infers. However, it dampens the variations that the Bayesian model infers: maxima are not as high in the SEIR model as in the Bayesian model, and minima not as low. These differences between the two approaches are significant and show that the SEIR model lacks in realism. In particular, it would seem adventurous to try and use it to obtain quantitative predictions about the future of the epidemic in France. However, we expect that despite these shortcomings, we can still use the SEIR model to obtain a qualitative picture of what the future of the epidemic may look like. 9

Migrations have a minor effect on the number of infections

In the following we examine the total number of infected individuals as of September 1 2020, just after the summer holidays, and at the end of the simulation, June 1 2021. Tables 1 and2 show that summer migrations have a very limited impact on the total number of infections. This limited impact at the national level is also observed at the regional level (data not shown). The impact seems slightly more important for low seasonality initially 1, but then this effect disappears at the end of the simulation 2. We also observed that migrations contribute to making the wave start slightly earlier (data not shown).

A second wave is observed in all parameter settings

The possibility of a second wave of infections depends on the strength of seasonality for SARS-Cov-2, and on its reproduction number R 0 . Fig. 5 shows the number of infectious individuals through time according to 2 different values for seasonality (low: 0.13, high: 0.5), and 3 different values for R average (1.1, 1.5, 2.0). Higher R average values result in earlier and steeper waves, because the epidemic spreads fast. When R average = 1.1, the wave is limited, peaking at 229 000 infectious individuals for low seasonality. Such a wave would be comparable to the first wave, which peaked at 307 000 individuals on March 29 according to [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF]'s model fitted on data up to June 29. Higher seasonality results in later and steeper waves. This is because, in our simulation, for a given R average , R 0 values of high seasonality scenarios are inferior to the R 0 values of low seasonality scenarios (Fig; 3). At the earliest, for low seasonality and high R average , the second wave peaks in October 2020; at the latest, for high seasonality and medium R average , the second wave peaks in January 2021 or even April 2021 for low R average .

Discussion

In this report we investigated the effect of several parameters on a model of the SARS-CoV-2 epidemic in France. We varied the value of the reproduction number, the value of the seasonality, and took into account summer holiday migrations in a SEIR model. Our choices of 1.1, 1.5, 2.0 for the future reproduction number of SARS-CoV-2 correspond to three hypotheses. In all cases, they assume that the sanitary measures currently in place, which include mandatory mask wearing in public transportation and at work, and the prohibition of large gatherings, are not sufficient to keep the reproduction number below 1.0, which means that the epidemic can spread. All choices are lower than the estimated value of the R average before the lockdown was enforced, which is around 3.0 [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF][START_REF] Salje | Estimating the burden of SARS-CoV-2 in France[END_REF][START_REF] Sofonea | Epidemiological monitoring and control perspectives: application of a parsimonious modelling framework to the COVID-19 dynamics in France. medRxiv[END_REF] , which conveys the idea that the sanitary measures do reduce the R average . The other parameters of the model were either obtained from estimates by [START_REF] Duchemin | Bayesian investigation of SARS-CoV-2-related mortality in France[END_REF])'s Bayesian model, or drawn from the literature. Variations in these parameters can affect the timing and the scale of the second epidemic wave (not shown), but unless extreme values are used, do not make it disappear.

We chose to use a deterministic SEIR model to quickly explore a range of scenarios. We do not attempt to take into account the uncertainty in most parameters, and do not attempt to model the inherent stochasticity of the epidemic process. The SEIR model assumes that all individuals behave the same way and have the same characteristics.

In particular, it assumes that all individuals in a given compartment and region have the same probability to be infected or to infect another individual. Models that allow heterogeneity between individuals have been proposed and can alter the predicted dynamics of the epidemic [START_REF] Gomes | Individual variation in susceptibility or exposure to SARS-CoV-2 lowers the herd immunity threshold[END_REF]. It is possible that such models may fit the dynamic of the epidemic better than our model, whose shortcomings are apparent Fig. 4. Future work may address the effect of migrations with such models. In the mean time we expect that our SEIR model can capture important aspects of the epidemic. The predictions of our model may provide insight into possible central trajectories of the epidemic, but cannot be used to infer e.g. a detailed timing of the progress of the epidemic.

Although the different French regions had different proportions of infectious individuals on July 1 2020 (Fig. 1), and the numbers of migrating individuals were large, we observed that holiday migrations have little effect on the epidemic (Tables 2,1). Our results are conservative, because we used an overestimate of the number of migrants during the summer holidays, which corresponds to the average number of migrants per year over 2014-2016, and not just during summer holidays. During migrations, at most 29 783 infectious individuals changed region on a given week (for R average = 2.0 and low seasonality). Compared to the 64 million individuals in France, these inputs were not sizeable enough to alter the dynamic of the epidemic in a meaningful way. This reasoning also explains why low seasonality results in a slightly larger impact of holiday migrations: in this condition higher summer R 0 values allow the epidemic to start an epidemic wave, meaning that the number of infectious migrants become less and less insignificant as time goes by. This shows that it is important that the number of infectious migrants remains low for migrations to have a limited impact on total numbers of deaths. However, small numbers of infectious migrants could still contribute to creating local clusters of infections. Investigating the likelihood of such clusters would require using an agent-based model with a realistic and detailed map of where individuals spend their holidays.

In a recent manuscript, [START_REF] Britton | Summer vacation and COVID-19: effects of metropolitan people going to summer provinces[END_REF] investigated the effect of summer migrations from a densely populated metropolitan area, Stockholm, to a popular summer vacation destination, the island of Gotland. They used a SIR model in which they introduced two types of visitors, long term visitors and short term visitors, assuming short term visitors have more contacts with the locals. They investigated a range of parameter values and found that migrations can increase the number of infections on the island. However, they also found that this increase depends upon the number of infectious migrants. Their results support our conclusion that low numbers of infectious migrants are key for migrations to have a limited impact on total numbers of deaths.

The reproduction number and the seasonality parameter have a large effect on the epidemic, and control the size, the spread and the date of a second wave. Even if the reproduction number is low at 1.1, a second wave is obtained. If the reproduction number is larger, the second wave happens earlier and is steeper and larger. This second wave is unsurprising given the low proportion of immunized individuals after the first wave, and given the historical precedent of the 1918 Influenza (Taubenberger & Morens, 13 2006).

Conclusion

We investigated the influence of summer holiday migrations and of various parametrizations for the seasonality and the reproduction number on the epidemic of SARS-CoV-2 in France. To this end, we modified a SEIR model to allow for migrations between regions during the summer holidays. We found that migrations of French individuals during summer holidays have almost no impact on the epidemic. This is a consequence of the lockdown which has been successful in decreasing the number of infectious migrants. Different seasonality and reproduction numbers above 1.0 all result in a second wave, whose peak date can vary between October 2020 and April 2021. Such a second wave would be avoided if the sanitary measures currently in place in France are sufficient to keep the reproduction number below 1.0. If the sanitary measures manage to keep the reproduction number low, for instance at 1.1 as in one of our simulations, the second wave is flattened and could be similar to the first wave. Such a flatter wave would be much more manageable for the French health system, hence it important to keep the reproduction number low.

Author contributions

LD provided code for running the Bayesian model to initialize the simulations. MP contributed to data analysis and acquisition. All authors discussed the ideas in the manuscript, and read and commented on the manuscript. BB directed the work, implemented the model and wrote the manuscript.

Figure 1 :

 1 Figure 1: Initial proportion of infectious individuals as of July 1 2020 according to Duchemin et al.'s Bayesian model. One unit µ corresponds to one infectious individual per 1 000 000.

Figure 2 :

 2 Figure 2: Migration balance during the summer holidays in continental France. The balance between outgoing and incoming migrants is color-coded between blue for regions with negative balance, and red for regions with positive balance.

Figure 3 :

 3 Figure 3: R t values combining seasonality and different values of R 0 . These different trajectories were used as input to the SEIR model.

  To validate our implementation of the SEIR model, we ran it between March 1 2020 and June 1 2021 and compared the number of infectious individuals through time to the numbers obtained by running (Duchemin et al. , 2020)'s Bayesian model over the same time interval. The SEIR model was initialized with the parameter values that the 8 Bayesian model has inferred for March 1 2020.

Figure 4 :

 4 Figure 4: Number of infectious individuals through time according to our SEIR model and according to a Bayesian model. Region-wise estimates from the Bayesian model are represented with a solid line on top of a shaded ribbon, and estimates from the SEIR model are represented with solid lines of matching colors.

Figure 5 :

 5 Figure 5: Number of infectious individuals through time. Scenarios with and without migration provide indistinguishable curves. 3 top lines: weak seasonality, R 0 at 1.1, 1.5, 2.0 from the top to the bottom. 3 bottom lines: strong seasonality, R 0 at 1.1, 1.5, 2.0 from the top to the bottom.

Table 1 :

 1 Total numbers of infections as of September 1 2020 for varying seasonality and R average values.

	10

Table 2 :

 2 Total numbers of infections as of June 1 2021 for varying seasonality and R average values.

Acknowledgements

The authors would like to thank (Neher et al. , 2020) for implementing their model in open source, which sped up our research, and Philippe Veber for fruitful discussions.