Tøke Camille Ameline
email: cameline8@gmail.com

Thomas Høye

ORCID Joseph James Bowden

Rikke Reisner Hansen

Oskar Liset

Pryds Hansen

Charlène Puzin

Philippe Vernon

Julien Pétillon

ORCID Toke Thomas Høye

Thomas Høye

Toke Thomas Høye

Elevational variation of body size and reproductive traits in high-latitude wolf spiders (Araneae Lycosidae)

Keywords: Arctic environment, biometry, altitude, environmental gradients, Pardosa, cocoon. Elevational Godthaabsfjord, Kob: Kobbefjord, Nar: Narsarsuaq, Far: Faroe Islands, Ice: Iceland. PE: postembryonic

Introduction

Life history traits are essential to determine the distribution of species and populations. In the Arctic, the relatively harsh environmental conditions can amplify biological phenomena witnessed in temperate regions, allowing us to observe changes in life history traits at smaller geographic scales [START_REF] Berteaux | Keeping pace with fast climate change: Can Arctic life count on Evolution?[END_REF][START_REF] Duarte | Abrupt climate change in the Arctic[END_REF]. For example, the same pattern of a given morphological trait might be visible along a shorter spatial gradient in the Arctic than in a temperate region. Because of this, the Arctic has become a major focus of ecological research in the context of climate change [START_REF] Høye | Arctic entomology in the 21st century[END_REF]. Furthermore, warming is much stronger in the Arctic and at high elevations [START_REF] Wang | Evidence of high-elevation amplification versus Arctic amplification[END_REF], and harsh environmental conditions generally lead to resource limitation [START_REF] Roff | Life history evolution[END_REF], causing allocation trade-offs: the Arctic is resource-poor compared to temperate regions (Bowden and Buddle 2012b;[START_REF] Ameline | Habitat specialization and climate affect arthropod fitness: a comparison of generalist vs. specialist spider species in Arctic and temperate biomes[END_REF]. Resource limitation is known to drive evolutionary strategies, such as the K, r and A strategies [START_REF] Greenslade | Adversity selection and the habitat template[END_REF][START_REF] Holm | Environmental restraints and life strategies: a habitat templet matrix[END_REF]. Individuals and populations are confronted with multiple environmental constraints along spatial gradients, like latitude or elevation. Abiotic factors such as temperature and growing season length vary across these gradients [START_REF] Conover | Countergradient variation in growth rate: compensation for length of the growing season among Atlantic silversides from different latitudes[END_REF], while biotic factors like predation risk also change [START_REF] Roslin | Higher predation risk for insect prey at low latitudes and elevations[END_REF]. Changes in elevation also trigger gradients in temperature lapse rate, short-wave radiation input, partial pressure of respiratory gases, precipitation (and related humidity), turbulence and wind speed (Hodkinson 2005[START_REF] Sundqvist | Community and ecosystem responses to elevational gradients: processes, mechanisms, and insights for global change[END_REF]). These changes can occur more abruptly across smaller spatial scales in elevation gradients compared to latitude gradients. Due to the rapid rate at which climate change is affecting northern ecosystems, it is essential to understand how populations respond to such climatic constraints [START_REF] Ohlberger | Climate warming and ectotherm body size -from individual physiology to community ecology[END_REF]. By combining latitudinal and elevational gradients, it may be possible to identify the environmental controls of life history trait variations.

Body size is a key trait related to the life history of individuals and can influence a broad range of factors from competition for mates to the amount of resources that an individual can allocate to reproduction. Body size variation along spatial gradients has been studied intensively, as body size is 'a continuously varying trait dependent on a variety of factors operating at different stages in an individual's life' [START_REF] Chown | Body size variation in insects: a macroecological perspective[END_REF]. The conditions under which an organism has developed are reflected in its body size. This makes body size a good labile trait for assessing phenotypic plasticity in response to climatic variability. Since the 19th century, many hypotheses were proposed as 'ecological rules' to try to explain observed variations of body size along geographical gradients. The most famous one is Bergmann's rule which proposes that the body size of endotherms increases with latitude (e.g. [START_REF] Atkinson | Why are organisms usually bigger in colder environments? Making sense of a life history puzzle[END_REF][START_REF] Blackburn | Geographic gradients in body size: a clarification of Bergmann's rule[END_REF][START_REF] Shelomi | Where are we now? Bergmann's rule sensu lato in insects[END_REF]. This rule states that a larger organism has a reduced surface to volume ratio, which minimizes heat loss. It has since been extended to include patterns in body size of ectotherms across latitudinal and elevational gradients [START_REF] Blanckenhorn | Bergmann and converse Bergmann latitudinal clines in arthropods: two ends of a continuum?[END_REF][START_REF] Chown | Body size variation in insects: a macroecological perspective[END_REF][START_REF] Shelomi | Where are we now? Bergmann's rule sensu lato in insects[END_REF], spiders: Entling et al. 2010). Bergmann's rule was set as a temperature effect while the converse Bergmann's rule has been set as a season length effect [START_REF] Sheridan | Shrinking body size as an ecological response to climate change[END_REF]. This extension to the rule supposes that shorter seasons at higher latitudes reduce the foraging period, resulting in a limited growth and development period [START_REF] Mousseau | Ectotherms follow the converse Bergmann's rule[END_REF]. A gradient of decreasing body size towards the poles should then be observable [START_REF] Blanckenhorn | Bergmann and converse Bergmann latitudinal clines in arthropods: two ends of a continuum?[END_REF]. Another rule, the countergradient variation is also linked to seasonal limitations and assumes that populations at higher latitudes would compensate for the reduced season length by growing faster [START_REF] Conover | Countergradient variation in growth rate: compensation for length of the growing season among Atlantic silversides from different latitudes[END_REF]. It is still a controversial subject and there is no clear explanation as to why body size increases or decreases when an organism is exposed to low temperatures [START_REF] Angilletta | Temperature, growth rate, and body size in ectotherms: fitting pieces of a life-history puzzle[END_REF]). Studies on ectotherms showed an increased body size of animals reared at lower temperatures, revealing a greater decrease of developmental rate than of growth rate (Forster et al. 2011a, b). This phenomenon is termed the 'developmental temperature-size rule' or TSR [START_REF] Atkinson | Temperature and organism size -a biological law for ectotherms?[END_REF]. There could also be a tradeoff between fecundity and a delayed maturity. Delayed maturation leads to a larger body size and an increase in fecundity [START_REF] Roff | Life history evolution[END_REF].

Wolf spiders (Lycosidae) are found in high densities in the Arctic (e.g. Bowden and Buddle 2012b). Among this group, female Pardosa carry their egg cocoons, or egg sacs, attached to their spinnerets, which allows for specimen-specific measurement of fecundity and quality of the offspring (e.g. [START_REF] Pétillon | Saline stress tolerance partly matches with habitat preference in ground-living wolf spiders[END_REF]. Clutch size has previously been considered a proxy for fecundity (e.g. [START_REF] Bowden | Fecundity and sexual size dimorphism of wolf spiders (Araneae: Lycosidae) along an elevational gradient in the Arctic[END_REF], as species of Pardosa in the Arctic are generally thought to produce only one clutch per lifetime and not to overwinter after reproduction (semelparous species, Bowden and Buddle 2012a). Consequently, we expect season length to play a larger role than temperature in the development of these spiders, compared to other species that reproduce multiple times in their lifetime or during multiple seasons. [START_REF] Høye | Climate change and sexual size dimorphism in an Arctic spider[END_REF] found in P. glacialis that body size varied mostly in relation to season length (timing of snowmelt).

Since clutch size is related to body size [START_REF] Pétillon | Plant invasion phenomenon enhances reproduction performance in an endangered spider[END_REF][START_REF] Puzin | Comparison of reproductive traits between two salt-marsh wolf spiders (Araneae, Lycosidae) under different habitat suitability conditions[END_REF]Hein et al. 2015[START_REF] Hein | Egg size versus egg number trade-off in the alpine-tundra wolf spider, Pardosa palustris (Araneae: Lycosidae)[END_REF], we also expect season length to have a greater impact on reproductive traits in Pardosa species. Egg volume commonly serves as a proxy for offspring quality, as bigger eggs result in bigger and fitter offspring (Hendrickx and Maelfait 2003).

In this study, we measure body size and reproductive traits along elevational gradients at six sites in Arctic and sub-Arctic areas. We aim to determine if closely related species (Pardosa spp.) display similar patterns in life history traits along these gradients. We test the following hypotheses: (1) As we expect season length to play a greater role than temperature in development, body size decreases with elevation, following a converse Bergmann's cline (spiders: latitude: Puzin et al. 2014, elevation: Høye and[START_REF] Bengtson | Natural History Museum Bern[END_REF][START_REF] Bowden | Fecundity and sexual size dimorphism of wolf spiders (Araneae: Lycosidae) along an elevational gradient in the Arctic[END_REF]Hein et al. 2015). This is the case in most terrestrial univoltine arthropods [START_REF] Horne | Temperature-size responses match latitudinal-size clines in arthropods, revealing critical differences between aquatic and terrestrial species[END_REF], and we expect similar constraints on the Pardosa species we studied. (2) As resources are more limited in high elevation, clutch size decreases with elevation; provided that female body size also decreases with elevation, as we know that reproductive outputs are linked to body size [START_REF] Bowden | Fecundity and sexual size dimorphism of wolf spiders (Araneae: Lycosidae) along an elevational gradient in the Arctic[END_REF]; and egg volume increases with elevation, revealing a quantity/quality trade-off in female fecundity in response to harsh environmental conditions (Hendrickx et al. 2003;[START_REF] Puzin | Comparison of reproductive traits between two salt-marsh wolf spiders (Araneae, Lycosidae) under different habitat suitability conditions[END_REF].

Material and Methods

Collected species

We collected four different species belonging to the genus Pardosa (Araneae, Lycosidae): P. hyperborea (Thorell 1872), P. palustris (Linnaeus 1758), P. glacialis (Thorell 1872) and P. furcifera (Thorell 1875), with contrasted habitat preferences, but most are thought to produce only one clutch per lifetime (semelparous species) (see Table 1 for life cycle information and Online Resource 1 for habitat data). Most Arctic and sub-Arctic Pardosa species have multiyear life cycles (at least two years), overwintering as young juveniles in their first winter and as sub-adults in their second winter [START_REF] Edgar | The life-cycle, abundance and seasonal movement of the wolf spider, Lycosa (Pardosa) lugubris, in Central Scotland[END_REF][START_REF] Pickavance | Life-cycles of four species of Pardosa (Araneae, Lycosidae) from the island of Newfoundland, Canada[END_REF][START_REF] Høye | Climate change and sexual size dimorphism in an Arctic spider[END_REF]Bowden and Buddle 2012b). Growing season spans from May-June to August-September depending on the species, location and sex (Table 1).

We collected P. hyperborea at three different sites in Greenland (Kobbefjord, Godthaabsfjord and Narsarsuaq), P. palustris in Iceland and in the Faroe islands, P. glacialis in Greenland (Zackenberg) and P. furcifera in Greenland (Kobbefjord) (Figure 1, Table 2). As we only collected three individuals of P. furcifera at high elevation, we did not test the effect of elevation on this species (Table 2). We hand collected over a small period of time (maximum of two weeks, Table 2) at each site. We collected thirty females with egg sacs at each point when possible (Table 2) within a relatively small area ('collection point'; < 50 m 2 ) and preserved them in 70% ethanol. We sampled either in 50 meter a.s.l. intervals along an elevational gradient or at low versus high elevations, depending on the elevational distribution of the species (Table 2, Online Resource 2). Lycosidae in boreal and Arctic regions are nonweb building species [START_REF] Marusik | The Greenland Entomofauna: An identification manual of insects, spiders and their allies[END_REF]. We believe that no migration was possible between collection points as they are separated by several kilometers, except for Kobbefjord (P. hyperborea) and Zackenberg (P. glacialis), where spiders were collected along a two kilometer regular slope, i.e. every 500 m (Online Resource 2).

Laboratory measurements

We estimated body size, fecundity and egg quality by prosoma width, clutch size, and mean egg volume, respectively in the laboratory with a EUROMEX trinocular stereomicroscope 45X, magnification 25X, connected with a CMEX5000 camera. Pictures and measurements were taken with the software EUROMEX ImageFocus v3.0. For each individual, we measured the following attributes: tibia length of leg I pair, length and width of the prosoma, clutch size (number of eggs) and mean length and width of ten randomly subsampled eggs (see Online Resource 3 for pictures). We calculated mean egg volume using the formula: Hendrickx and Maelfait (2003). We divided embryonic stages into six categories: I: only cells are visible; II: legs are visible; III: the prosoma is developed; IV: the prosoma is hatched and no measurement is possible; V: the postembryonic individual is hatched within the egg sac; VI: the postembryonic individual has emerged from the egg sac (Online Resources 3 and 4). We used this last stage when the collection in ethanol made the postembryonic individuals (1 st instar larvae) exit the egg sac.

𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉𝑉 = 𝜋𝜋 6 × 𝑉𝑉𝑉𝑉𝑙𝑙𝑙𝑙𝑙𝑙ℎ × 𝑤𝑤𝑤𝑤𝑤𝑤𝑙𝑙ℎ 2 following
Since we always collected a female with an intact egg sac, we assume not to have lost any individuals. A large proportion of the egg sacs (34 % of all individuals) contained eggs, which had hatched when dissected (stages IV to VI). This bias is however difficult to avoid given that sampling cannot be done with exact synchrony in the field. Due to the variation in the development of eggs, we considered mean egg volume as a less reliable trait than clutch size, and refrained from using a global variable connecting clutch size and mean egg volume. We did all measurements on all four species, with a total of 775 females. Of the 775 individuals in total, 14 cocoons (13-low, 1-high elevation) of P. hyperborea were found to be parasitized from the Narsarsuaq (South Greenland) site. No parasitism was found in the other three species. This represented an incidence of 1.8% parasitism for all data combined, but 18% of the population of P. hyperborea collected from the low site at Narsarsuaq. Teneral adults were identified as Gelis sp. (Hymenoptera, Ichneumonidae, Cryptinae) (Online Resource 3, Figures 2.1 to 2.3).

Analyses

In order to get a representative measurement of body size, correlations were assessed between all of our size measurements. Prosoma width was correlated with prosoma length, left tibia length and right tibia length by 99% (Pearson correlation test, t = 170.0, df = 763, p < 0.0001), 95% (t = 83.8, df = 738, p < 0.0001) and 95% (t = 84.9, df = 732, p < 0.0001), respectively. Hence, we chose to use prosoma width as a proxy of body size (as used in [START_REF] Jakob | Estimating fitness: a comparison of body condition indices[END_REF]. Henceforth we refer to prosoma width as body size if not stated otherwise. First, the effects of elevation on body size and on reproductive traits were assessed. Sampling differences (Table 2) necessitated that we tested for the effects of elevation on species and sites, separately. Second, the effect of body size on reproductive traits was assessed for our entire data set and for each species separately. All measurements were normally distributed (Online Resource 5). Linear models were applied, LLSR (Linear Least Squares Regression) in the case of continuous explanatory variables, ANCOVA (Analysis of Covariance) in the case of both continuous and categorical explanatory variables and factors. When testing for the effect of elevation, we included collection point as a random effect to avoid site-specific effects using linear mixed models (function 'lmer' in R package 'lme4', [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF].

For testing the effects of elevation on reproductive traits, we always include female body size as a covariate. Significance of explanatory variables and factors was tested via an ANOVA test using type III sums of squares (package 'car': Fox and Weisberg 2011), i.e., without taking the order of the explanatory variables/factors into account. Finally, normal distribution of the residuals of the model was tested through a distribution histogram and a quantilequantile plot. All models are presented in Tables 3 and4. All means are presented with standard error (package 'RVAideMemoire': Hervé 2015) and all statistics were conducted using the R environment for statistical computing v. 3.3.2 (R Development Core Team 2016).

Results

The effect of elevation on body size and reproductive traits Mean body size (prosoma width) of Pardosa hyperborea, P. palustris, P. glacialis and P. furcifera were 1.85 ± 0.006 mm (n = 465), 2.17 ± 0.013 mm (n = 90), 2.60 ± 0.012 mm (n = 151) and 3.05 ± 0.020 mm (n = 59), respectively. Mean body size of P. hyperborea was larger in Kobbefjord with mean sizes of 1.95 ± 0.0096 mm (n = 178), 1.80 ± 0.0071 mm (n = 166) in Godthaabsfjord and 1.78 ± 0.0077 mm (n = 123) in Narsarsuaq, (ANOVA, p < 0.0001, TukeyHSD test, p < 0.0001). Elevation had no effect on the body size of P. glacialis (p = 0.56) or of P. hyperborea in Godthaabsfjord (p = 0.72) or Kobbefjord (p = 0.44, Table 3, Figure 2a). There was a positive effect of elevation on body size for P. hyperborea in Narsarsuaq (p = 0.049, Table 3, Figure 2a). At the lower elevation, mean body size was 1.77 ± 0.009 mm (n = 52) compared to 1.80 ± 0.013 mm (n = 71) at the higher elevation. Similarly, elevation had a positive effect on P. palustris body size (p = 0.0092, Table 3, Figure 2a): mean body size was 2.09 ± 0.019 mm (n = 30) at the lower elevation compared to 2.23 ± 0.023 mm (n = 30) at the higher elevation.

Elevation did not have an effect on clutch size for P. glacialis (p = 0.40), P. hyperborea in Godthaabsfjord (p = 0.19) and Narsarsuaq (p = 0.065) and for P. palustris (p = 0.12) (Table 3, Figure 2b). For P. hyperborea in Kobbefjord, elevation did have a negative effect on clutch size (p = 0.00039, Table 3, Figure 2b): clutch size at the lowest elevation was 39 ± 1 (n = 60) compared to 27 ± 1 (n = 30) at the highest elevation. Elevation did not have any effect on mean egg volume for all species and sites: P. glacialis (p = 0.46), P. hyperborea in Godthaabsfjord (p = 0.39), Kobberfjord (p = 0.73) nor P. palustris (p = 0.96) (Table 3, Figure 2c). This effect was not tested for P. hyperborea in Narsarsuaq as only 5 females carried unhatched eggs.

The effect of body size on reproductive traits

There was an effect of species (p < 0.0001), a positive effect of female body size (p < 0.0001) and an effect of their interaction on clutch size (p < 0.0001) (Table 4), meaning that the effect differed between species. Specifically, the positive effect of female body size on clutch size was stronger for larger species (Figure 3, regression slopes: P. hyperborea: 33.04 ± 2.26, P. palustris: 47.53 ± 7.28, P. glacialis: 82.93 ± 7.20, P. furcifera: 81.02 ± 12.77). There was no effect of species (p = 0.33), body size (p = 0.15) nor their interaction (p = 0.44) on mean egg volume (Table 4), meaning that mean egg volume varied independently from female body size.

Discussion

Body size patterns along elevational gradients among sites and species

For Pardosa palustris and P. hyperborea at one collection site, female body size increased between low and high elevation, however female size in P. glacialis and P. hyperborea at the two other sites did not respond to elevation. Hence, it appears that body size of P. palustris and P. hyperborea are more sensitive to elevation than P. glacialis and that the smaller species are more sensitive to elevational gradients. This result contradicts our first hypothesis, which was that spiders would be smaller at higher elevations. Similarly, Hein et al. (2015) did not find any general body size pattern along elevational gradients. Conversely, Høye and Hammel (2010) did find a decrease in body size on female P. glacialis at high elevation.

Body size patterns along elevation could be species-specific (Høye and Hammel 2010;[START_REF] Bowden | Fecundity and sexual size dimorphism of wolf spiders (Araneae: Lycosidae) along an elevational gradient in the Arctic[END_REF] or may depend on the body size of the species in question [START_REF] Blanckenhorn | Bergmann and converse Bergmann latitudinal clines in arthropods: two ends of a continuum?[END_REF]. Such variation could arise from the combined effect of life-history traits, habitat, latitude and movement as we discuss below.

We observed different patterns among species that can be triggered by differences in lifehistory traits (see Table 1). We know that P. glacialis individuals take multiple years to reach adult size [START_REF] Høye | Climate change and sexual size dimorphism in an Arctic spider[END_REF][START_REF] Marusik | The Greenland Entomofauna: An identification manual of insects, spiders and their allies[END_REF]: Table 1), but smaller species may require less time to reach adulthood. We assume that the larger species P. furcifera also has a longer lifecycle, although this is not documented. These species could thus increase their development time to be able to grow larger, trading delayed maturity for a larger body size, in spite of the possible disadvantages of a longer development time. This suggests that larger Arctic Pardosa species might have developed a distinct strategy to tackle environmental challenges.

Ecological trade-offs have also been described in temperate forest wolf spiders, e.g. risk of freezing vs. ability to grow at low temperatures [START_REF] Whitney | Tradeoff in two winter-active wolf spiders: increased mortality for increased growth[END_REF]). However, phylogeography [START_REF] Gür | Why do Anatolian ground squirrels exhibit a Bergmannian size pattern? A phylogenetic comparative analysis of geographic variation in body size[END_REF], adaptation [START_REF] Stillwell | Are latitudinal clines in body size adaptive?[END_REF] or biological interactions (e.g., [START_REF] Callaway | Positive interactions among alpine plants increase with stress[END_REF] might also help explain size responses to elevation.

Temperature is known to decrease approximately 0.5-0.55°C with an increase of 100 m in elevation [START_REF] Marshall | Near-surface-temperature lapse rates on the Prince of Wales Icefield, Ellesmere Island, Canada: implications for regional downscaling of temperature[END_REF], while growing season becomes shorter with elevation [START_REF] O'neill | Warm' Tundra: Atmospheric and near-surface ground temperature inversions across an alpine treeline in continuous permafrost, western Arctic, Canada[END_REF]. We know that climate change is stronger in the Arctic and at high elevations [START_REF] Wang | Evidence of high-elevation amplification versus Arctic amplification[END_REF], which causes earlier snowmelt and thus relatively longer growing seasons [START_REF] Semmens | Early snowmelt events: detection, distribution, and significance in a major sub-arctic watershed[END_REF]. It is estimated that a 1°C increase in mean annual temperature would result in a range change of ~167 m in elevation and ~145 km in latitude [START_REF] Jump | The altitude-for-latitude disparity in the range retractions of woody species[END_REF]. It has been shown in a large variety of terrestrial arthropods that larger species are less affected by warming than smaller species [START_REF] Forster | Warming-induced reductions in body size are greater in aquatic than terrestrial species[END_REF]. This could allow large Pardosa species to become even larger in response to climate change.

Intraspecifically, it remains difficult to predict the effects of increasing temperatures along elevational gradients, as elevation-body size clines often appear to deviate from temperaturebody size cline measured in the laboratory [START_REF] Horne | Insect temperature-body size trends common to laboratory, latitudinal and seasonal gradients are not found across altitudes[END_REF].

Local habitat variation could also drive the patterns in body size along elevational gradients, as habitats change rapidly over elevational gradients. Indeed, body size trends in similar groups have been shown to be affected by local habitat variation [START_REF] Høye | Climate change and sexual size dimorphism in an Arctic spider[END_REF][START_REF] Bowden | Habitat-specific effects of climate change on a low-mobility Arctic spider species[END_REF]. Previous studies of Arctic spider assemblages have also suggested that habitat is a key determinant in species turnover across elevation in these ecosystems (Bowden andBuddle 2010, wet vs. dry habitats: Ernst et al. 2016;Hansen et al. 2016b;[START_REF] Cameron | Seasonal change and microhabitat association of Arctic spider assemblages (Arachnida: Araneae) on Victoria Island (Nunavut, Canada)[END_REF][START_REF] Høye | Elevation modulates how Arctic arthropod communities are structured along local environmental gradients[END_REF].

In our study, we could not disentangle the effect of species from that of latitude because all species were not present at all sampling sites, and this is why we ran separate models for each species and sites. As a consequence, latitudinal effects could mask or supress elevational effects. We found distinct patterns of body size and reproductive traits along separate elevational gradients in the same species (P. hyperborea) sampled at nearby sites (the Kobbefjord population was sampled between 40 km and 100 km away from the Godthaabsfjord populations during the same period in the same year (Table 2), which suggests that latitude is not the primary factor that affects elevational patterns in our study species. We suggest that variation in microhabitat influenced our results, while inter-annual climatic variation could also have played a role. We found our largest species, P. furcifera, in South Greenland, although only three individuals were collected at high elevations. Further sampling of this species or other large species such as P. groenlandica are necessary to test the hypothesis that larger species are less sensitive to environmental changes, regardless of latitude. We sampled P. glacialis in Zackenberg (North-East Greenland), although it is also found at lower latitudes in South-West Greenland but at low abundances and only at high elevations (Hansen et al. 2016a). Adaptation to harsher climatic conditions at higher latitudes could also explain why we do not observe body size nor reproductive traits responses to elevation in P. glacialis. Some Pardosa species have been shown to be able to actively migrate to maximum distances of 500 m per day (P. pullata, P. prativaga, [START_REF] Richter | Differences in breeding and motility between Pardosa pullata (Clerk) and Pardosa prativaga (L. Koch), (Lycosidae, Araneae) in relation to habitat[END_REF], 50 m in 10 min (P. lapidicina, [START_REF] Morse | Orientation and movement of wolf spiders Pardosa lapidicina (Araneae, Lycosidae) in the intertidal zone[END_REF], 280 m per day (P. monticola in moss dunes, [START_REF] Bonte | Patch quality and connectivity influence spatial dynamics in a dune wolfspider[END_REF], or 40-150 m over several months (P. pullata, P. chelata, Hallander 1967). We know that wolf spiders are able to migrate in response to environmental cues [START_REF] Kraus | Seasonal habitat shift in an intertidal wolf spider: proximal cues associated with migration and substrate preference[END_REF], and such shifts to suitable microhabitats may constitute a protective behaviour [START_REF] Schaefer | Winter ecology of spiders (Araneida)[END_REF]). Since we do observe a pattern of body size in Kobbefjord in P. hyperborea, we suppose that migration does not play a significant role there. P. glacialis was collected along an elevational gradient of just 200 meters as well, but we do not observe a pattern of body size, so perhaps this environmental gradient was not sufficient to elicit an effect in this species.

Interspecific differences in reproductive traits responses to elevation

We found that clutch size and egg volume varied independently of elevation for P. glacialis and P. palustris while for P. hyperborea in Kobbefjord, clutch size decreased with elevation, but egg volume did not show any pattern. Similarly, [START_REF] Hein | Egg size versus egg number trade-off in the alpine-tundra wolf spider, Pardosa palustris (Araneae: Lycosidae)[END_REF] found no variation of the eggs size/number trade-off in P. palustris along an elevational gradient. Our results show that P. hyperborea responded more strongly to elevation than the other two species. This supports our hypothesis of increased sensitivity of smaller Pardosa species to environmental changes, since P. hyperborea was the smallest species we collected. The absence of a trend between egg volume and body size can also indicate that the volume or mass of individual eggs is influenced more by a combination of body condition and environmental factors than clutch size (e.g. Bowden and Buddle 2012b). As egg volume does not vary along elevational gradients, it appears unnecessary for the eggs to be larger in colder environments.

Interspecific variation in the relationship between body size and clutch size

For all the species we studied, clutch size increased with body size, while egg volume varied independently of body size. This shows that intraspecifically, larger individuals could increase their fecundity more efficiently than smaller individuals. It has been shown among arthropods that clutch size generally increases with body size (Fox and Czesak 2000), including wolf spiders [START_REF] Pétillon | Plant invasion phenomenon enhances reproduction performance in an endangered spider[END_REF][START_REF] Puzin | Comparison of reproductive traits between two salt-marsh wolf spiders (Araneae, Lycosidae) under different habitat suitability conditions[END_REF]Hein et al. 2015[START_REF] Hein | Egg size versus egg number trade-off in the alpine-tundra wolf spider, Pardosa palustris (Araneae: Lycosidae)[END_REF]. This emphasizes that a larger body size is beneficial for females because it enables them to increase their fecundity and ultimately their fitness.

From an interspecific point of view, the larger the species, the higher the rate at which clutch size increases with body size. Hence, larger species could increase their fecundity more efficiently than smaller species. According to the Metabolic Theory of Ecology (MTE), having a larger body size allows one to retain more energy because it is processed slower [START_REF] Brown | Toward a metabolic theory of ecology[END_REF]). This energy can thus be invested in a higher, albeit usually delayed reproductive output compared to smaller species. Within species, this theory can also explain why larger individuals can invest more in their reproductive output. Here, the reproductive output is associated with clutch size while egg volume is associated with the quality of offspring.

Our observation that egg volume was not well predicted by body size in wolf spiders mirrors findings by others [START_REF] Brown | Clutch size and offspring size in the wolf spider Pirata sedentarius (Araneae, Lycosidae)[END_REF]Hendrickx and Maelfait 2003;Hein et al. 2015). Indeed, investment in individual egg mass or volume does appear to be better explained by body condition than body size in this group (e.g., Bowden and Buddle 2012b).

Conclusion

We set out to assess the impact of environmental variation, employed through the use of elevational gradients, on body size and reproductive traits in Arctic wolf spiders. Due to the relative harshness of the environmental conditions in this region, we expected that elevational gradients would serve as local seasonal and temperature gradients. We found that smaller species at the respective sites showed more sensitivity to environmental changes. Second, intra-and interspecifically, larger individuals could increase their fecundity more efficiently than smaller individuals/species. Together, these results suggest that large body size is advantageous in the Arctic environment for these species. Third, it appears that females do not produce larger eggs at higher elevations. This could be explained by the fact that the high variability of alpine microclimatic conditions makes it possible for spiders to find suitable conditions even at higher elevations. Our results also show that there are no large-scale general patterns in species-specific responses to environmental changes, and that microclimatic conditions along with inter-annual climatic variation probably play a crucial role. We emphasize the need for further studies to measure, along with large-scale gradients, local environmental variables such as soil moisture, temperature, shrub dominance, snow-cover and growing season length. [i] highest activity: mid-June to late-July at high latitudes [d] males and females: June to August [f] one clutch per lifetime [a] egg sacs: July and August [f,i] , hatching:

August [i] one or two winters before maturation, females are thought not to overwinter after reproduction [a,b] overwinters in soil cracks and under stones

[c] P. hyperborea 2-year life span, inferred from related species [b,g] males: late May to early August, females: June to September [f] July to (at least) September [b] highest activity: mid-June to late-July at high latitudes [d] one clutch per lifetime [a] August to September [b] July [e] egg sacs observed in July and August [f] one or two winters before maturation, females are thought not to overwinter after reproduction [a,b] P. furcifera 2-3 year life span, inferred from related species [b,g] * males: mid-June to July, females: mid-June to August

[f] highest activity: mid-June to late-July at high latitudes [d] one clutch per lifetime [a] egg sacs: July, hatched: mid-August [f] one or two winters before maturation, females are thought not to overwinter after reproduction [a,b] P. palustris 2-3 year life span, [b,g,h,j] May and June [f] 1-3 clutches per lifetime [a,j] egg sacs: June [f] one or two winters before maturation, females are thought not to overwinter after reproduction [a,b,h] 551

Fig. 1

 1 Fig.1 Map of the collecting sites of Pardosa species. From North to South: Zac: Zackenberg (74.467N, 20.566W), Ice: Iceland (64.700N, 21.500W), God: Godthaabsfjord (64.144N, 51.474W), Kob: Kobbefjord (64.117N, 51.350W), Far: Faroe Islands (62.041N, 6.877W), Nar: Narsarsuaq (61.183N, 45.371W). Latitude and longitude in Decimal Degrees WGS84.

Fig. 2

 2 Fig.2 Effect of elevation on body size, clutch size and mean egg volume. a. body size, b. clutch size, c. mean egg volume. Sample size was too low for P. hyperborea in Narsarsuaq to apply a statistical test (n = 5). Significant values are written in bold. Applied models are detailed in Table3

Fig. 3

 3 Fig.3 Linear regression results of clutch size on female size. Formulas are written as follows: 𝑦𝑦 = 𝑎𝑎 + 𝑏𝑏 * 𝑥𝑥. Pardosa hyperborea: LLSR, clutch = -31.68 + 33.04*size, F 445 = 212.7, R 2 = 0.32, p < 0.0001, P. palustris: LSSR, clutch = -60.52 + 47.53*size, F 88 = 42.67, R 2 = 0.32, p < 0.0001, P. glacialis: LLSR, clutch = -141.42 + 82.93*size, F 149 = 132.7, R 2 = 0.47, p < 0.0001, P. furcifera: LLSR, Clutch = -159.28 + 81.02*prosoma width, F 57 = 30.49, R 2 = 0.34, p < 0.0001

Table 1

 1 Life cycle information currently known for the collected species. a: Bowden and Buddle 545 2012a; b: Pickavance 2001; c: Marusik 2015; d: Bowden and Buddle 2010; e: Hammel and Nickel 546 2008; f: Dondale and Redner 1990; g: Buddle 2000; h: inferred from related species in similar latitudes 547 e.g. Edgar 1971; i: Høye et al. 2009; j: Hein et al. 2018. * We propose that P. furcifera could have a 548 similar life-span to P. glacialis, as they present similar sizes

	549
	550
	Species	Life cycle	Development season	Reproduction Hatching period	Overwintering
	P. glacialis	2-3 years

Acknowledgments

This work was part of the SPACEWOLF project (Spatial gradients in physiological adaptation and life history traits of Arctic wolf spiders) led by Julien Pétillon and Philippe

556

* pools that were not used for models assessing elevation's effects on body size and reproductive traits.

557

v pools that were not used for the model assessing elevation's effects on mean egg volume due to variation in development

558

stages. Thorell 1872 Holarctic [20] , Nearctic, Arctic, Greenland (G) [21] not specified 0-860 [13,14,15] P. hyperborea Thorell 1872 Holarctic [20] , low Arcticboreal [21] , Greenland (S, SW north to Disko, SE north to 70°50'N) [21] rare psychrophilous [7] 100-1360 [3,5,12,13,14,15,18] P. furcifera Thorell 1875 Canada, Alaska, Greenland (S, SW [21] ), Iceland [20] : trans-Nearctic arcto-alpine range [10, 21] not specified 0-1360 [13,14,15] P. palustris Linnaeus 1758

560

Holarctic [20] generalist [7] (widespread, especially in cold climates) [1] up to 2500 [2,7,9,11] Species Found habitats P. glacialis in heaths [21] and sedge meadows, among lichens and low shrubs, and on beaches or bare soil on Arctic tundra [6] , transition boreal forest to Arctic tundra [14] , sub-Arctic pine forest area, birch woods, bogs [21] and treeless fells [15] P. hyperborea open forests and moors [19] , birch wood [2,15] , sphagnum bogs (in boreal zones [8] ), low-lying spruce forests, rocky hillsides, among lichens in Arctic and alpine tundra and on needle mats in forests ofjack pine, black spruce, and balsam fir [6] , transition boreal forest to Arctic tundra [14] , sub-Arctic pine forest area, birch woods, bogs and treeless fells [15] , south-facing dry slopes [21] , heath habitats [22] P. furcifera spruce-pine forests, in willow clumps near the timberline, and on moss and stones in alpine tundra [6] , sub-Arctic pine forest area, birch woods, bogs and treeless fells [15] , transition boreal forest to Arctic tundra [14] , on rather moist ground and in Salix herbacea snow beds as well as in the moss layer of heaths adjoining brooks. Also in herb fields on slopes and in Salix glauca thickets [21] P. palustris open, dry areas [19,17,1] , tidal drifts on a saltwater beach, along grassy margins of streams, above timberline in mountain passes [6] , rich grassland, poor grassland, ash woodland [9] , tundra, mire [16,18] ,tall herb meadow, grass meadow, birch wood (open ground and sparse woods) [1,2,3] , pastures, heaths [1] , bogs [8] , braided rivers [3] , grass and dwarf shrub heaths and open infield meadows [4] G: global, S: south, W: west, N: north, E: east

Polar Biology

Camille Ameline (corresponding author), Toke Thomas Høye, Joseph James Bowden, Rikke Reisner Hansen, Oskar Liset Pryds Hansen, Charlène Puzin, Philippe Vernon, Julien Pétillon

Corresponding author: cameline8@gmail.com Online Resource 3 Fig. 1.1 to 1.6: body measurements and species' prosomas. Fig. 2.1 to 2.3: parasitoid developmental stages. Fig. 3.1 to 3.12: eggs developmental stages. Arrows represent the measurements. We used a EUROMEX trinocular stereomicroscope 45X connected with a CMEX5000 camera. Pictures and measurements were taken with the software Euromex ImageFocus v3.0.

mm


Vineeth S Varma

Irinel-Constantin Morȃrescu
email: constantin.morarescu@univ-lorraine.fr

Dragan Nešić

Open multi-agent systems with discrete states and stochastic interactions

Keywords: Multi-agent systems, Markov chains, open systems, consensus

We consider a multi-agent system in which agents arrive and depart from a network randomly as a Bernoulli process. Each agent that is active in the network must decide between two actions represented by 0 or 1. Each active agent then observes the action of a random neighbour and updates its preference towards a certain action. New agents that arrive into the network are activated with a random preference and action. This means that the notion of consensus in the standard sense can no longer be applied and instead, we provide conditions under which majority action preservation occurs when the number of agents is arbitrarily large. This property will imply that a large fraction of the active agent population will retain their action almost surely.

I. INTRODUCTION

The analysis of multi-agent systems received an increasing attention during the last decades. They have a wide range of applications covering robotics [START_REF] Bullo | Distributed Control of Robotic Networks[END_REF], [START_REF] Mesbahi | Graph theoretic methods in multiagent networks[END_REF], power networks [START_REF] Chow | Time scale modeling of sparse dynamic networks[END_REF], [START_REF] Bıyık | Area aggregation and time-scale modeling for sparse nonlinear networks[END_REF] and opinion dynamics [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF], [START_REF] Friedkin | Social influence and opinions[END_REF]. In the framework of multi-agent systems, the most studied problem is the consensus achievement [START_REF] Jadbabaie | Coordination of groups of mobile autonomous agents using nearest neighbor rules[END_REF], [START_REF] Hendrickx | Convergence of type-symmetric and cut-balanced consensus seeking systems[END_REF], [START_REF] Moreau | Stability of multiagent systems with time-dependent communication links[END_REF], which captures various coordination behaviors. Some of the most cited features leading to the usage of multi-agent systems are their flexibility and scalability. While the multi-agent systems literature mainly treats the case of fixed dimension networks, the size of physical networks can increase or decrease in time. Indeed, a typical example is the one of a social network whose size is time-varying and often very large but it does not increase unbounded. On the other hand, technological networks may also suffer from failure/repairing of some agents leading to modifications of the active subsystems in the network. This motivates us to analyze here the behavior of open multi-agent systems, in which agents keep arriving and/or leaving the network during the evolution in time.

Open multi-agent systems have been empirically studied in [START_REF] Török | Opinions, conflicts, and consensus: modeling social dynamics in a collaborative environment[END_REF], [START_REF] Iniguez | Modeling social dynamics in a collaborative environment[END_REF] in the framework of opinion dynamics. In computer science, there exist specific architectures coping with the possibility of agents joining and leaving the system [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF]. A formal mathematical analysis of open multi-agent systems with deterministic arrivals and departures has been recently provided in [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF], in which, all the dynamics are deterministic, and agents are only allowed to arrive, and do not depart.

As for [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF], [START_REF] Varma | Modeling stochastic dynamics of agents with multi-leveled opinions and binary actions[END_REF], a precise motivation of the present work is given by opinion dynamics in social networks. Instead of imposing pre-determined or periodic time instants for arrivals like [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF], we assume that arrivals and departures occur with certain probabilities. Another important difference is that we consider random gossiping over the time-varying complete graph of active agents, with each agent having access only to a quantized version of the state of random neighbours, and not to the exact value. An important idea that we introduce in this paper, which is also highly relevant for opinion dynamics, is the one of almost sure consensus. This also represents a major difference with respect to [START_REF] Varma | Modeling stochastic dynamics of agents with multi-leveled opinions and binary actions[END_REF] in which we study a closed system (i.e., the number of agents is fixed) that always leads to consensus. We note that gossiping under quantized communication for closed systems has been studied in [START_REF] Frasca | Average consensus on networks with quantized communication[END_REF]. Beside the fact that we consider open systems, in our model, the state of any agent can evolve only inside a discrete/finite set of values which seems more appropriate for many opinion dynamics applications.

The main challenges that we face in this work are the timevarying dimension of the system and the absence of classical "convergence". To overcome the first challenge, we consider that there exists a fixed, but arbitrarily large dimension of the system denoted by N . Consequently, at any instant of time t, only a fraction of the N agents are active, and their number characterizes the real size of the network, while the others are inactive and artificially added to facilitate the analysis. We assign a precise state value to active agents that form the network at time t while inactive agents have a generic inactive state referred to as NA.

The second challenge is that consensus in the classical sense is impossible to achieve due to repeated arrivals with a different state, and this is indeed one reason why not many works have theoretically studied consensus in open systems. Therefore, we propose and analyze a different notion of consensus, which requires that a large fraction of the active population maintains a common action (i.e., the quantized version of the state) almost surely. The main objective of this work is to provide conditions under which, the proposed dynamics stabilizes the consensus points almost surely. Our results are suitable for application to several scenarios where consensus of discrete actions are required such as synchronizing the wireless channel with which mobile devices communicate, opinion dynamics, etc.

The rest of the paper is organized as follows. Section II introduces the main concepts and formulates the problem under consideration. Section II-B mathematically formulates the system dynamics as a Markov chain. The main results concerning the equilibria of large-scale multi-agent systems under the proposed stochastic dynamics are presented in Section III. The results of our numerical studies and simulations are illustrated in Section IV. The paper ends with some concluding remarks and perspectives.

Notation. We use R for the set of real numbers, R + the set of non-negative reals and N = {0, 1, . . . } the set of non-negative integers. The indicator function is given by 1 S (x) which takes the value 1 when x ∈ S and 0 otherwise. Throughout the paper the function x represents the integer part of x + 1 2 . The sign function which is 1 when x > 0, 0 if x = 0 and -1 if x < 0 is denoted by sgn(x). We also use O(•) for the big O notation describing the limiting behavior. As our consensus dynamics has stochastic elements, we also introduce Pr(A) to denote the probability of event A, and E[•] for the expectation of a random variable. The notation Pr(A|B) is used for the the probability of event A conditioned by the event B. We also use almost surely to describe events that happen with probability 1 and define the corresponding stability notion for any stochastic process p(t).

Definition 1 (Almost surely stable): If there exists * > 0 s.t.,

Pr |p(t + 1) -p * | < |p(t) -p * A | < = 1 (1)
for all < * A , then we say that p * is almost surely stable for the stochastic dynamics of p(t).

II. PROBLEM STATEMENT A. System dynamics

We consider a discrete-time multi-agent system with t ∈ N denoting the time instants. We use V := {1, 2, . . . , N } to denote the set of agents. At any given time instant t, only a subset V A (t) ⊆ V of agents are active in the network. We use V I (t) := V \ V A (t) to denote the set of inactive agents. Let us also introduce a fixed value L ∈ 2N + 1 and the set L = {l 0 , l 1 , . . . , l L } with l k = k L , ∀k ∈ {0, 1, . . . , L}. When an agent is inactive, we assign the state NA, i.e. for all i ∈ V I (t), x i (t) = NA. We call the extended state space

L E = L ∪ {NA}.
We assign to any agent i ∈ V A (t) a hidden state/preference x i (t) ∈ L and a visible action based on its state as

q i (t) = x i (t) . (2)
It is worth noting that, when L = 1, the state and the action sets are identical. Since we consider an open system, any agent which is active at time t i.e., i ∈ V A (t), may become inactive at t + 1 with a probability δ ∈ (0, 1). Similarly, any inactive agent i ∈ V I (t) may activate at t + 1 with a probability γ. This is:

Pr x i (t + 1) = NA | x i (t) ∈ L = δ, Pr x i (t + 1) ∈ L | x i (t) = NA = γ,
On top of this, at any time t such that |V A (t)| ≥ 2 (i.e., more than one agent is active), an agent i ∈ V A (t) may have only one neighbour with index J i (t) which is randomly chosen among all the other active agents:

Pr(J i (t) = j) = 1 |V A (t)| -1 , ∀i ∈ V A (t), j ∈ V A (t) \ {i}.
(3) Moreover, the agent i ∈ V A (t) will update with probability β ∈ (0, 1 -δ) its state x i (t) by observing the action of its neighbour J i (t) ∈ V A (t) \ {i} as follows:

Pr x i (t + 1) = φ(x i (t), q Ji(t) (t)) x i (t) ∈ L = β (4) where φ( , q) = + 1 L sgn (q -) (5)
for any ∈ L, q ∈ {0, 1}. This model is inspired by [START_REF] Martins | Discrete opinion models as a limit case of the coda model[END_REF] in which a Bayessian update rule was proposed in the context of opinion dynamics. Additionally, classical discrete-state multi-agent dynamics like the Ising model, which can be seen as a special case of (4) when L = 1, do not lead to consensus in open systems as we show later in Remark 2.

In other words, the system dynamics is defined by random gossiping which indicates that any active agent can become a neighbour of i with identical probabilities. Finally, the active agent may chose to retain its preference if the above events don't occur, i.e.,

Pr x i (t + 1) = x i (t) | x i (t) ∈ L = 1 -δ -β, Pr x i (t + 1) = x i (t) | x i (t) = NA = 1 -γ. (6)
When an inactive agent enters the network and becomes active, it must decide on a certain state. We consider that agents chose a random state with identical probabilities, i.e.,

Pr x i (t + 1) = x i (t) = NA = γ L + 1 (7)
for any ∈ L. Note that any agent which becomes inactive does not retain any memory of its preference while it was active. We can use ( 2)-( 7) to describe the proposed algorithm for achieving majority consensus in action/decision. When |V A (t)| ≥ 1, define p -(t) as the fraction of the active population with action 0, and p + (t) as the fraction with action 1, calculated as

p + (t) = n∈V A (t) q n (t) |V A (t)| (8)
and p -(t) = 1 -p + (t). Throughout the paper we also use p A (t) to denote the fraction of active agents, i.e.,

p A (t) = |V A (t)| N = N n=1 1 L (x n (t)) N (9)
Next, at any t such that p A (t) > 0, we define

p (t) = N n=1 1 { } (x n (t)) N n=1 1 L (x n (t)) (10)
for any ∈ L, which denotes the fraction of active agents with state . Note that we can also write

p + (t) = n∈V A (t) q n (t) |V A (t)| = L k= L+1 2L p l k (t) (11)

B. Markov model

In order to give a complete mathematical model of the open system that takes into account all the stochastic phenomena associated with the activation/deactivation of agents as well as possibility of update or not, we formulate the dynamics as a Markov process in which ∀i ∈ V, ∈ L:

                                 Pr(x i (t + 1) = + sgn(1-) L | x i (t) = ∈ L) = β Pr( x Ji(t) (t) = 1) = β j∈V A (t)\{i} qj (t) |V A (t)|-1 , Pr(x i (t + 1) = + sgn(-) L | x i (t) = ∈ L) = β Pr( x Ji(t) (t) = 0) = β j∈V A (t)\{i} 1-qj (t) |V A (t)|-1 , Pr(x i (t + 1) = NA | x i (t) ∈ L) = δ, Pr(x i (t + 1) = x i (t) | x i (t) ∈ L) = 1 -β -δ, Pr(x i (t + 1) = | x i (t) = NA) = γ L+1 , Pr(x i (t + 1) = NA | x i (t) = NA) = 1 -γ. ( 12
)
Remark 1: In the framework of opinion dynamics, the first two equations in [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF] express that any individual in the social network (active agent) updates its state according to the update rule (4) with probability β. The third equation formalize the fact that individuals of the network can leave it at time t (i.e., any active agent at time t becomes inactive at time t+1) with probability δ. The fourth equation describe the possibility of one individual in the network to not interact and preserve its opinion at time t (i.e., an active agent preserves its state) with probability 1 -δ -β. The last two equations express the possibility of individuals outside network to join it (choosing the state according to (7)) or stay inactive, respectively.

We can see from (12) that x(t) is a Markov process as the transition probabilities only depend on its current state. However, note that x(t) ∈ L N E . Therefore, even when L = 1, the Markov chain has 3 N states, and 5 N states when L = 3. In order to find the stationary state distribution of such a process would involve finding the eigenvalues of an (L + 2) N × (L + 2) N matrix which is infeasible for large N (even N > 10). Indeed, many modern networks have a huge number of agents in practice. These two reasons motivate us to analyze the system in the large scale limit, i.e., when N → ∞.

The main objective of this work is to analyze the behavior of dynamics [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF] and show that a majority of active agents in the network will reach a consensus in their action almost surely. That is, by implementing the proposed algorithm, we obtain that ∃ * > 0 s.t. ∀ ∈ (0, * ),

         Pr p + (t + 1) < p + (t) < = 1 or Pr p -(t + 1) < p -(t) < = 1 (13)

III. LARGE SCALE NETWORK ANALYSIS

For convenience, we define p := (p l0 , p l1 , . . . , p l L ) T as the population vector. This allows us to define the following.

Definition 2 (Equilibrium in population): We say that (p * A , p * ) is an equilibrium in population for dynamics (12) if ∃t ≥ 0 such that (p A (t), p(t)) = (p * A , p * ) for all t > t . Note that the p * A denotes the population fraction of active agents and p * describes the distribution of the active population among the states in L with ∈L p = 1. Therefore both p * A and p * are required to describe the configuration of the system in terms of population. Furthermore, (p * A , p * ) is almost surely stable for the overall system if p * A and p * are respectively almost surely stable for the dynamics p A (t) and p(t).

Proposition 1: When N → ∞, we have

p A (t + 1) = (1 -δ)p A (t) + γ(1 -p A (t)) (14)
almost surely, which results in an almost surely globally exponentially stable (GES) equilibrium at p * A = γ γ+δ . Proof: Since the activation and deactivation probabilities are given by γ and δ respectively, and these probabilities are independent of the actions/preferences, we can always write

Pr(x n (t + 1) ∈ L) = (1 -δ)1 L (x n (t)) + γ1 {NA} (x n (t))
(15) for any n ∈ V based on [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF], resulting in

E[ N n=1 1 L (x n (t + 1))] = (1 -δ) N n=1 1 L (x n (t)) +γ N n=1 1 {NA} (x n (t)) = N ((1 -δ)p A (t) + γ(1 -p A (t))) (16
) Note that the random variable 1 L (x n (t + 1)) has a variance greater than 0 and less than 1. Since N → ∞, we can use the central limit theorem [START_REF] Billingsley | Convergence of probability measures[END_REF] to conclude that the sampled average p A (t+1) will converge almost surely to the expectation, i.e., (1 -δ)p A (t) + γ(1 -p A (t)). The dynamics ( 14) can be rewritten as

p A (t + 1) -p A (t) = -(γ + δ)p A (t) (17)
where p A (t) := p A (t) -γ γ+δ . Since γ + δ ∈ (0, 2), the origin is GES for the deterministic dynamics (17) implying γ γ+δ is GES for [START_REF] Varma | Modeling stochastic dynamics of agents with multi-leveled opinions and binary actions[END_REF]. However, as [START_REF] Varma | Modeling stochastic dynamics of agents with multi-leveled opinions and binary actions[END_REF] only occurs almost surely, the equilibrium γ γ+δ is almost surely stable for the p A dynamics.

This result allows us to characterize the population of active agents while at equilibrium. Since this equilibrium is independent of the fraction of agents with specific preference or action states, any equilibrium of the system will have p * A = γ γ+δ . We characterize the dynamics of p , for all ∈ L under this equilibrium as follows.

In the following, we consider that γ and δ depend on N . We recall that N is fixed but arbitrarily large and our analysis is based on the fact that large scale networks behavior is approximated by N → ∞. Therefore, we will add supplementary assumptions on the behavior of γ and δ when N → ∞.

Proposition 2: If N → ∞ and δ = O(γ), then when p A (t) = p * A = γ γ+δ (equilibrium of active agent population), we have

p l0 (t + 1) = (1 -δ)p l0 (t) -βp l0 (t)p + (t) + δ L+1 +βp l1 (t)p -(t) p l k (t + 1) = (1 -β -δ)p l k (t) + βp l k-1 (t)p + (t) + δ L+1 +βp l k+1 (t)p -(t) p l L (t + 1) = (1 -δ)p l L (t) -βp l L (t)p -(t) + δ L+1 +βp l L-1 (t)p + (t) (18)
for all k ∈ {1, . . . , L -1} almost surely.

Proof: Let us suppose that the number of active agents is at equilibrium i.e., p A (t) = γ γ+δ . As we assume that δ = O(γ) one has that lim

N →∞ j∈V A (t)\{i} q j (t) |V A (t)| -1 = p + (t) and
we can write

lim N →∞ Pr(x n (t + 1) = l k ) = (1 -δ -β)1 l k (x n (t)) +β1 l k-1 (x n (t))p + (t) + γ L+1 1 NA (x n (t)) +β1 l k+1 (x n (t))p -(t) (19
) for any k ∈ {1, . . . , L -1} based on [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF]. Using the central limit theorem as done before in Proposition 1, we can evaluate

N n=1 1 l k (x n (t)) N = (1 -β -δ)p l k (t)p A (t) + γ L+1 (1 -p A (t)) + βp l k-1 (t)p + (t)p A (t) +βp l k+1 (t)p -(t)p A (t) (20)
almost surely. Dividing by p A (t) > 0 on both sides, we have the left hand side becoming p l k (t + 1) and using γ 1-p A (t) p A (t) = δ one obtains [START_REF] Ising | Contribution to the theory of ferromagnetism[END_REF] for all ∈ {l 1 , . . . , l L-1 }. We can similarly evaluate for p l0 (t + 1) and p l L (t + 1), but by keeping in mind that φ(l 0 , 0) = l 0 and φ(l L , 1) = l L .

Proposition 2 converts the stochastic dynamics of our Markovian system [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF] into a deterministic population dynamics (18) that will occur with probability 1 when N → ∞. Next, we characterize the equilibrium points of (18) when δ = O(1/N ) and γ = O(1/N ).

Theorem 1:

When N → ∞, δ = O(1/N ), γ = O(1/N ), N γ
γ+δ → ∞, L ≥ 3 and p A (t) = p * A , the system (12) has exactly two almost surely locally exponentially stable equilibrium points at p * S1 = (1, 0, . . . , 0) T and p * S2 = (0, . . . , 0, 1) T . The system also allows for one unstable equilibrium point at

p * U = 1 L + 1 • (1, 1, . . . 1) T .
Proof: See Appendix A. From Theorem 1, we can see that taking L ≥ 3 will result in almost all active agents with a common action at equilibrium when δ = O(1/N ) and γ = O(1/N ). Since Theorem 1 provides local exponential stability of (18) around p * S1 , we can write

∃ * > 0 s.t. ∀ < * , Pr (p (t + 1) < |p (t) < ) = 1 (21)
for all ∈ L \ {l 0 } which shows that the first condition of objective ( 13) is satisfied. Considering the other locally stable point p * S2 , we can also prove ( 13) is satisfied for the other condition.

Although these results are derived for N → ∞, we perform simulations and show that that the probability of retaining consensus is close to 1 when N is large but finite.

Remark 2 (Case of L = 1): The case of L = 1 which is the only case excluded from Theorem 1 (recall L is odd) implies that x n (t) = q n (t) as L = {0, 1} and this results in a dynamics similar to an Ising model used in [START_REF] Ising | Contribution to the theory of ferromagnetism[END_REF]. We observe that equation ( 18) can be simplified into

p l0 (t + 1) = (1 -δ)p l0 (t) -βp l0 (t)p l L (t) + βp l L (t)p l0 (t) +δ/2 = (1 -δ)p l0 (t) + δ/2 p l L (t + 1) = (1 -δ)p l L (t) -βp l0 (t)p l L (t) + βp l L (t)p l0 (t) +δ/2 = (1 -δ)p l L (t) + δ/2
(22) which is almost surely GES at p * U and which results in p + (t) = p -(t) = 0.5. Therefore taking L = 1 ensures that there will be no consensus at equilibrium.

IV. NUMERICAL RESULTS

In this part, we numerically illustrate the theoretical results presented in the previous Section. These simulations show that as far as the arrivals/departures occur sufficiently rarely with respect to the dimension of the network (N δ and N γ are small) we obtain majority action consensus.

For our simulations we take L = 3, i.e. L = {0, 1/3, 2/3, 1} and we fix β = 0.5. In our first set of simulations we plot the population of agents with action 1, i.e., p + (t) vs t with δ = γ = 0.01 for various values of N . In Figure 1a, when N = 10, we see how consensus is sometimes perturbed by the activation of new agents. When N sufficiently large as in Figure 1b, consensus is maintained for a very long duration of over 10 5 time instants.

In Figure 2, we plot the population with action 1 against time for N = 1000 and various values of δ = γ. When δ = 0.01, the simulation results are close to the theoretical behavior predicted for N → ∞ and δ = O(1/N ). When δ = 0.2, the influence of the agents who get activated or deactivated is dominant resulting in random actions, and there is no majority action preservation. Interestingly, even with δ = 0.1, about 75% of the agents maintain a common action over a long duration.

V. CONCLUSION

We introduced and analyzed an open multi-agent system with discrete states and quantized information. Basically, the system is split into two subsets: active agents and inactive ones. Each active agent is characterized by a state belonging to a finite set while the inactive agents are just artificially introduced to preserve the dimension of the overall system. The system dynamics is modeled as a gossiping in which the active agents access only the action (quantization of the state) of neighbours. At any time the network can be affected by arrivals and departures of some agents (i.e., some active agents can become inactive and vice-versa). These events occur randomly as a Bernoulli process. We show that, in this framework, classical consensus can be replaced by a weaker formulation: a large fraction of the active population has almost surely reached a common action.

  (a) Dynamics of each agent when N = 10. Agent activation is marked by a circle and deactivation with a square. At around t = 560, the majority agreement value at 0 is perturbed to become 1.(b) Population dynamics, i.e. fraction of agents with action 1 vs time when N = 100 and the almost surely stable equilibrium (defined when N → ∞) persists for a very long duration (much more than 10 5 steps).

Fig. 1 :

 1 Fig. 1: Agent dynamics for δ = 0.01.

Fig. 2 :

 2 Fig.2: Plotting p + (t) vs t for various values of δ. We observe that even when δ = 0.1, i.e. about 100 agents are being activated or deactivated in the network at each time instant, a majority of agents hold the same action for a long period of time.

This work was partially supported by the Australian Research Council under the Discovery Project SP170104099 and the French CNRS PEPS project YPSOC.

APPENDIX

A. Proof of Theorem 1

Proof: In Proposition 2 we have shown that the population dynamics associated with the Markovian system [START_REF] Carrascosa | Multiagent Systems, "Service oriented mas: an open architecture[END_REF] can be given by [START_REF] Ising | Contribution to the theory of ferromagnetism[END_REF] almost surely. Since we have δ = O(1/N ) and γ = O(1/N ), as N → ∞, δ, γ → 0, (18) will be transformed into p l0 (t + 1) = p l0 (t) -βp l0 (t)p + (t) + βp l1 (t)p -(t) p l k (t + 1) = (1 -β)p l k (t) + βp l k-1 (t)p + (t) +βp l k+1 (t)p -(t) For the stability analysis of (23), we perform a linear analysis of the dynamics around these equilibria by studying the Jacobian matrix. Denote by g k (p) the dynamics of the population with preference state l k ∈ L, i.e. p l k (t + 1) = g k (p). For the rest of the proof, for ease of exposition, we omit the argument t from p (t).

If we denote the Jacobian elements by J i,j , where J i,j = ∂gi ∂p l j

, then for all 0 < i ≤ L-1 2 , and for all L+1 2 ≤ j < L, we have:

We also have ,

For all i, j ∈ {1, 2, . . . , L -1} such that |i -j| > 1 we have

when j ≤ (L -1)/2 and

when j ≥ (L + 1)/2. Next, we have

For all i ∈ {1, 2, . . . , L -1}, we have

where k = i+1 if i+1 ≤ (L-1)/2 and k = i-1 otherwise; and

where k = i+1 if i-1 ≤ (L-1)/2 and k = i-1 otherwise. We can study the stability of an equilibrium point by looking at the eigenvalues of Jacobian matrix evaluated at the equilibrium points, i.e.,

which is a triangular matrix. This matrix therefore has the eigenvalues 1 and 1 -β. This can be easily verified by evaluating the determinant of the J(p * S1 ) -λI. The eigenvector corresponding to the eigenvalue of 1 is in fact (1, 0, . . . ). Since we are working on the unit simplex (as ∈L p = 1), the equilibrium p * S1 is locally exponentially stable within the unit simplex. By symmetry, we can also show the same for the equilibrium p * S2 . However, the dynamics (23) occurs almost surely and is not deterministic, which results in these two equilibria being almost surely stable. Now, we study the Jacobian at the other equilibrium p * U . Here, p + = p -= 0.5 and the J at this point can be evaluated. Its first column (J 0,0 , J 1,0 , . . . , J L,0 ) T is given by

The columns j for j ∈ {1, . . . , (L -1)/2} are of the form

where 1 + β L+1 + β 2 is the diagonal term of the Jacobian. The j-th column for j ∈ {(L + 1)/2, . . . , L -1} is of the form

) T Finally, the L-th column is given by

The above matrix is such that each column has exactly one element which is -β L+1 either at the first row (after column index is more than (L -1)/2) or at the last row. We can evaluate the sum of each column as

when L ≥ 3 for all j ∈ {0, 1, . . . , L}. We can show that one of the eigenvalues of the matrix J(p * U ) is 1 + β L+1 (L -2). To do this, we look at

Using the property that adding a scalar times a row to another row does not change the determinant, we replace the first row with the sum of all rows, and this results in the first row becoming all zeroes since we know from (24) that the sum of each column is 0 by adding -1-(L-2) β L+1 once to each row. Since J(p * U ) has an eigenvalue 1 + (L -2) β L+1 > 1 when L ≥ 3, we see that p * U is unstable when L ≥ 3.