Claire Jacquet
email: claire.jacquet@wanadoo.fr

Luca Carraro

Florian Altermatt

Meta-ecosystem dynamics drive the spatial distribution of functional groups in river networks

Keywords: Dendritic network, macroinvertebrates, metaecosystem, resource transportation, river continuum concept, riverine ecosystem, OCN

The meta-ecosystem concept provides a theoretical framework to study the effect of local and regional flows of resources on ecosystem dynamics. Meta-ecosystem theory has hitherto been applied to highly abstract landscapes, and meta-ecosystem dynamics in real-world landscapes remain largely unexplored. River networks constitute a prime example of meta-ecosystems, being characterized by directional resource flows from upstream to downstream communities and between the terrestrial and the aquatic realm. These flows have been thoroughly described by the River Continuum Concept (RCC), a seminal concept in freshwater ecology, stating that observed spatial variations in the relative abundances of invertebrate functional groups reflect systematic shifts in types and locations of food resources, which are in turn determined by the physical attributes of river reaches. Hence, the RCC represents a solid conceptual basis for determining how changes in landscape structure and resource flows will translate into local and regional changes in community composition. Here, we develop and analyse a riverine meta-ecosystem model inspired by the RCC, which builds upon a physically-based landscape model of dendritic river networks. We show that the spatial distributions and regional biomass of invertebrate functional groups observed in stream communities are determined by the spatial structure and scaling attributes of dendritic river networks, as well as by specific rates of resource flows. Neglecting any of these aspects in modelling river meta-ecosystems would result in different community patterns. Moreover, we observed that high rates of resource flow, for example due to river anthropization, have a negative effect on the regional biomass of all functional groups studied, and can lead to cascading extinctions at the meta-ecosystem scale. Our work paves the way for the development of physically-based meta-ecosystem models to understand the structure and functioning of real-world ecosystems.

Introduction

Ecosystems are open to flows of materials, organisms, and energy, and understanding the effects of these flows on ecosystem structure and functioning is a central goal of ecological research [START_REF] Polis | Food webs at the landscape level[END_REF]. By integrating the local production and movement of resources into metacommunity theory, the meta-ecosystem framework allows investigating feedback processes between community and resource dynamics across spatial scales [START_REF] Loreau | Meta-ecosystems: A theoretical framework for a spatial ecosystem ecology[END_REF][START_REF] Massol | Linking community and ecosystem dynamics through spatial ecology[END_REF]Gounand et al., 2018a;[START_REF] Guichard | Meta-Ecosystem Dynamics[END_REF]. Seminal work on meta-ecosystems has been developed on simplified landscapes [START_REF] Loreau | Meta-ecosystems: A theoretical framework for a spatial ecosystem ecology[END_REF][START_REF] Loreau | Spatial flows and the regulation of ecosystems[END_REF], often represented by two-patch ecosystems (Gravel et al., 2010a,b;[START_REF] Marleau | Nutrient flows between ecosystems can destabilize simple food chains[END_REF]. This simple representation allowed an analytical investigation of the interplay between local and regional flows of matter and their implications for community dynamics [START_REF] Loreau | Meta-ecosystems: A theoretical framework for a spatial ecosystem ecology[END_REF]Gravel et al., 2010a;[START_REF] Massol | Linking community and ecosystem dynamics through spatial ecology[END_REF][START_REF] Leroux | Dynamics of Reciprocal Pulsed Subsidies in Local and Meta-Ecosystems[END_REF]. More recent studies focusing on larger spatial networks outlined the importance of spatial structure and movement rates of organisms and materials to promote meta-ecosystem stability [START_REF] Marleau | Meta-ecosystem dynamics and functioning on finite spatial networks[END_REF][START_REF] Gravel | Stability and complexity in model metaecosystems[END_REF].

However, one enduring limitation of current meta-ecosystem models is the abstract representation of the landscape (Gounand et al., 2018a), which is often described by random or Cartesian spatial networks [START_REF] Marleau | Meta-ecosystem dynamics and functioning on finite spatial networks[END_REF][START_REF] Gravel | Stability and complexity in model metaecosystems[END_REF]. In contrast, the physical structure of real-world landscapes constrains organisms' movement and resource flows and is likely to influence the spatial distributions of resources and the organisms that feed upon them [START_REF] Leroux | Subsidy hypothesis and strength of trophic cascades across ecosystems[END_REF]Harvey et al., 2017a;[START_REF] Schmitz | Animals and the zoogeochemistry of the carbon cycle[END_REF][START_REF] Montagano | The strength of ecological subsidies across ecosystems: a latitudinal gradient of direct and indirect impacts on food webs[END_REF][START_REF] Harvey | Metaecosystem dynamics drive community composition in experimental, multi-layered spatial networks[END_REF].

River networks constitute a prime example of meta-ecosystems, as documented by a large body of literature in freshwater ecology assessing the critical influence of directional resource flows, both from the surrounding terrestrial environment and from upstream river reaches, on the composition of local stream communities (e.g., [START_REF] Bartels | Reciprocal subsidies between freshwater and terrestrial ecosystems structure consumer resource dynamics[END_REF]; [START_REF] Soininen | Toward more integrated ecosystem research in aquatic and terrestrial environments[END_REF]; [START_REF] Abelho | Litter movement pathways across terrestrial-aquatic ecosystem boundaries affect litter colonization and decomposition in streams[END_REF]). Evidence that composition of local riverine communities is linked with the river ecosystem in its entirety was prominently pointed out in the River Continuum Concept (RCC) [START_REF] Vannote | The river continuum concept[END_REF], a cornerstone concept in freshwater ecology. The RCC proposes that commonly observed spatial variations in the relative abundances of major functional groups of organisms along a longitudinal river gradient reflect systematic shifts in types and locations of food resources, which are in turn determined by the hydrological attributes of river reaches (e.g., stream width, depth and velocity) and by the energetic constraints of upstream communities.

The organismal groups considered are freshwater invertebrates that are generally clustered in five functional groups (grazers, shredders, collectors, filterers and predators), which are of high relevance with respect to the biodiversity and functioning of riverine ecosystems [START_REF] Anderson | Detritus processing by macroinvertebrates in stream ecosystems[END_REF][START_REF] Wallace | Stream invertebrate productivity linked to forest subsidies: 37 stream-years of reference and experimental data[END_REF][START_REF] Harvey | Regulation of the functional structure of aquatic communities across spatial scales in a major river network[END_REF]. The RCC specifically predicts that shredders should be the most abundant functional group in small river reaches, which are strongly shaded by the surrounding vegetation and receive a large input of dead organic matter from falling leaves. The abundance of grazers is expected to peak at mid-sized streams, where light penetration into the stream is highest, stimulating the development of primary producers.

Light penetration, and therefore primary production, is limited by water depth and turbidity in large rivers, reducing the abundance of grazers. Finally, filterers and collectors should be the most abundant groups in large rivers where the most abundant resources is fine particulate organic matter, which is a by-product of leaf consumption by shredders and is delivered into the water column from upstream communities.

The RCC is likely the most influencing conceptual framework in freshwater and stream ecology and is among the most commonly cited works in this field (totalling >6000 citations in the Scopus database), prompting empirical and theoretical works assessing its strengths and limitations (see [START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF] for a recent review). Although the parallel between meta-ecosystem theory and the mechanisms formulated in the RCC has been acknowledged in several recent studies [START_REF] Massol | Linking community and ecosystem dynamics through spatial ecology[END_REF]Gounand et al., 2018a;[START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF][START_REF] Harvey | Metaecosystem dynamics drive community composition in experimental, multi-layered spatial networks[END_REF], a formal integration of the RCC within a meta-ecosystem model is still lacking. Next to the explicit predictions mentioned above, the RCC further implicitly predicts that changes in the hydrological attributes of a river network should translate into local and regional changes in the composition of invertebrate communities. However, this prediction has never been investigated formally and it remains unclear how changes in resource input and transport will affect the spatial distribution of invertebrate functional groups in a river network. Furthermore, the RCC essentially describes community patterns along the main stem of a river, hence disregarding the contribution of river dendritic structure on the spatial patterns observed in river meta-ecosystems [START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF].

In this study, we develop a spatially explicit meta-ecosystem model for river systems inspired by the RCC to investigate the effect of meta-ecosystem dynamics on the functional composition of stream communities. We specifically address the following research questions: (i) What is the contribution of the dendritic structure of river networks on the spatial distribution of functional groups described in the RCC? (ii) How do changes in resource flow rate influence the spatial distribution and regional biomass of functional groups in river networks? To do so, we make use of a physically-based model of dendritic river networks, which expresses how stream characteristics (e.g., water discharge, stream width) vary across a river system as a function of drainage area (i.e., the portion of land over which precipitation is drained towards a given river cross-section).

We compare the spatial distributions of functional groups between a complex dendritic river network (see [START_REF] Carraro | Optimal Channel Networks accurately model ecologicallyrelevant geomorphological features of branching river networks[END_REF] and three alternative landscapes based on different assumptions with respect to dendritic structure and scaling of hydrological variables.

We demonstrate that the spatial patterns described in the RCC can only emerge from a metaecosystem model that accounts for the dendritic structure of river networks, as well as for changes of its hydrological attributes in the downstream direction. We then analyse different scenarios of resource spatial dynamics and show that increased rates of resource flows have a negative impact on the regional biomass of all the functional groups studied and can lead to extinctions at the meta-ecosystem scale.

Methods

Meta-ecosystem structure and dynamics

We considered a riverine meta-ecosystem composed of a set of local ecosystems, each of them defined by a river reach and its surrounding terrestrial area, which are spatially connected via resource flow across the river network. A river reach corresponds to an uninterrupted stretch of river in which abiotic conditions can be assumed as constant, and hence constitutes the fundamental unit of a river network. In our model, each local ecosystem is composed of four non-living resources, a group of primary producers and five groups of consumers, corresponding to the trophic groups and resource types commonly found in freshwater ecosystems [START_REF] Vannote | The river continuum concept[END_REF][START_REF] Wallace | The role of macroinvertebrates in stream ecosystem function[END_REF][START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF][START_REF] Larsen | Testing the River Continuum Concept with geostatistical stream-network models[END_REF].

First, grazers G feed on primary producers P (algae or aquatic rooted vascular plants), which are supported by nutrients N (e.g., nitrogen) and light (Fig. 1). Second, three independent functional groups feed on different forms of particulate organic matter derived from decomposed terrestrial leaf litter: shredders S, collectors C, and filterers F feed on coarse (CPOM, > 1 mm in size), fine (FPOM, 50 µm-1 mm), and ultra-fine (UPOM, < 50 µm) particulate organic matter, respectively. Finally, predators R feed on other groups of consumers, that is G, S, C and F (Fig. 1).

The resources present in a local ecosystem originate both from local terrestrial inputs of N, CPOM and UPOM and from upstream reaches via hydrological transport of N, CPOM, FPOM and UPOM. Consequently, local ecosystems are connected through hydrological transport of resources, and the composition of functional groups in upstream ecosystems determines not only the local resource consumption, but also the amount of resources available for the ecosystems situated downstream. Note that [START_REF] Vannote | The river continuum concept[END_REF] originally included both filterers and collectors in the same functional group (i.e., collectors), while a distinction between collectors (feeding on FPOM) and filterers (feeding on UPOM) has been operated in more recent literature (see e.g. [START_REF] Larsen | Testing the River Continuum Concept with geostatistical stream-network models[END_REF][START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF]. Hence, the densities of collectors described in [START_REF] Vannote | The river continuum concept[END_REF] should be compared with the sum of collector and filterer densities in our model.

We focused our analysis on an equilibrium state for the meta-ecosystem, and did not consider the temporal fluctuation of resource inputs, such as seasonality in resource availability or stream flow. We further hypothesized that living organisms only move within river reaches and do not disperse across local ecosystems. Our meta-ecosystem model is expressed by the following set of ordinary differential equations:

dR i dt = R α R R i (G i + S i + C i + F i ) -µ R R i -β R R 2 i ; (1a) dG i dt = G α G G i P i -α R R i G i -µ G G i -β G G 2 i ; (1b) dS i dt = S α S S i CPOM i -α R R i S i -µ S S i -β S S 2 i ; (1c) dC i dt = C α C C i FPOM i -α R R i C i -µ C C i -β C C 2 i ; (1d) dF i dt = F α F F i UPOM i -α R R i F i -µ F F i -β F F 2 i ; (1e) dP i dt = P α P l i P i N i -α G G i P i -µ P P i -β P P 2 i ; (1f) dN i dt = δ N V i n ∑ j=1 w ji N j Q j -N i Q i + φ N,i V i -α P l i P i N i -λ N N i ; (1g) dCPOM i dt = δ CPOM V i n ∑ j=1 w ji CPOM j Q j -CPOM i Q i + φ CPOM,i V i -α S S i CPOM i -λ CPOM CPOM i ; (1h) dFPOM i dt = δ FPOM V i n ∑ j=1 w ji FPOM j Q j -FPOM i Q i + (1 -S )α S S i CPOM i -α C C i FPOM i + -λ FPOM FPOM i ; (1i) dUPOM i dt = δ UPOM V i n ∑ j=1 w ji UPOM j Q j -UPOM i Q i + φ UPOM,i V i -α F F i UPOM i -λ UPOM UPOM i , ( 1j
)
where subscript i identifies a river reach. System Eq. ( 1) is based on a Lotka-Volterra formulation of trophic interactions (including self regulation terms, following [START_REF] Barbier | Pyramids and cascades: a synthesis of food chain functioning and stability[END_REF]), and on mass balance equations (Eqs. ( 1g)-( 1j)) for the dynamics of resources in river networks. For the living compartments X = {R, G, S, C, F, P} (Eqs. ( 1a)-( 1f)), the rate of change in biomass density X depends on biomass gain owing to feeding (assuming direct dependence of feeding rate on biomass production) and biomass loss due to metabolism, predation and intra-group competition for resources (density dependence). The parameters related to the living compartments X are feeding rate α X , efficiency of resource assimilation X , strength of intra-group competition β X , and mortality rate µ X .

To include the effect of light availability on nutrient uptake of producers, we expressed the realized assimilation rate of producers as α P l i , where l i is a site-specific light limitation factor derived from physical principles and by following the downwelling irradiance concept [START_REF] Fasham | A nitrogen-based model of plankton dynamics in the oceanic mixed layer[END_REF][START_REF] Davies-Colley | Predicting light penetration into river waters[END_REF]. To derive l i , we assumed that the irradiance of photosynthetically active radiation above the canopy is constant across the river system, while spatial variations in l i are solely determined by variations in river geometry (i.e., width and depth-see Supporting Information for details). Where available, parameters were chosen in agreement with established evidence on feeding behaviour for the various functional groups.

Detailed information on parameter values is available in Table S1.

For the resource compartments Y = {N, CPOM, FPOM, UPOM} (Eqs. Eqs. ( 1g)-( 1j)), the rate of change of resource density Y depends on mass gain from local terrestrial inputs and inputs from upstream reaches via hydrological transport; and mass loss due to consumption by living organisms, downstream hydrological transport, and resource degradation and deposition.

The parameters related to resources Y are the flux of local terrestrial inputs φ Y,i , the relative downstream velocity of resource Y with respect to water δ Y (i.e., if δ Y = 0.5, resource Y travels downstream half as fast as water), and the rate of resource loss due to processes other than consumption λ Y (i.e., degradation or deposition). Q i and V i represent the water discharge and the water volume of reach i, respectively; w ji is the generic entry of the adjacency matrix (w ji = 1 if river reach j drains into i and 0 otherwise); n is the total number of reaches constituting the river network.

Following [START_REF] Vannote | The river continuum concept[END_REF]; [START_REF] Marks | Revisiting the fates of dead leaves that fall into streams[END_REF], we assumed that FPOM is a by-product of CPOM consumption by shredders (Fig. 1), therefore, as a first approximation, inputs of FPOM originate from the aquatic environment only. Importantly, Eqs. (1g-j) outline how the hydrological attributes of river networks, that is water discharge Q i and volume V i , influence the density of resources available locally for the functional groups. In particular, both water discharge and volume determine the velocity at which resources are transported across a given reach, while water volume also influences the concentration of local terrestrial inputs.

Resource spatial dynamics in river networks

We generated a large, virtual river network (so-called optimal channel network, OCN) with the R-package OCNet [START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF][START_REF] Carraro | OCNet: Optimal Channel Networks[END_REF] in order to study changes in the abundance density of functional groups along a gradient of physical conditions specific to river networks. OCNs are structures that reproduce the topological connectivity and scaling features of real river networks [START_REF] Rinaldo | Evolution and selection of river networks: Statics, dynamics, and complexity[END_REF] and are well suited to study riverine ecological processes [START_REF] Carraro | Generation and application of river network analogues for use in ecology and evolution[END_REF][START_REF] Carraro | Optimal Channel Networks accurately model ecologicallyrelevant geomorphological features of branching river networks[END_REF]. We built the OCN (Fig. 2a) on a square lattice spanning an area of 5625 km 2 partitioned into 3688 reaches in order to obtain a gradient of physical conditions sufficiently large to reproduce the predictions of the RCC, which apply to river systems spanning a wide range of stream sizes. For instance, this roughly corresponds to the size of river Rhone before entering Lake Geneva (5238 km 2 ) or river Rhine before entering Lake Constance (6299 km 2 ).

By using a threshold area of 1 km 2 , the resulting total river length is 3854 km (see details in the Supporting Information).

We made several assumptions regarding the rate of hydrological transport of resources in order to incorporate the verbal arguments of the RCC into our meta-ecosystem model. In particular, we assumed CPOM to be transported downstream at a low rate (δ CPOM = 0.01) compared to water because of the large size of its constituents, which likely induces clogging [START_REF] Vannote | The river continuum concept[END_REF][START_REF] Wallace | The role of macroinvertebrates in stream ecosystem function[END_REF]. Conversely, we assumed FPOM to travel downstream at an intermediate rate (δ FPOM = 0.5), while UPOM and N are transported at a high rate (i.e., same velocity as water: δ UPOM = 1 and δ N = 1) [START_REF] Cushing | Transport dynamics of fine particulate organic matter in two Idaho streams[END_REF][START_REF] Wallace | The role of macroinvertebrates in stream ecosystem function[END_REF].

We derived water discharge Q i and water volume V i (Fig. 2c) across all river reaches of the network based on drainage area values (Fig. 2b) via the scaling relationships of [START_REF] Leopold | The hydraulic geometry of stream channels and some physiographic implications[END_REF]. Drainage area corresponds to the portion of land over which precipitation is drained towards a given location (Fig. 2a). As a universal geomorphological feature, drainage area is the master variable controlling the physical and hydrological characteristics of a given reach [START_REF] Leopold | Fluvial Processes in Geomorphology[END_REF][START_REF] Rodriguez-Iturbe | Fractal River Basins. Chance and self-organization[END_REF]. We therefore used drainage area to describe the positioning of a reach within the river network and illustrated how water volume (Fig. 2c), concentration of terrestrial inputs (Fig. 2d), and light availability (Fig. 2e) change along a gradient of drainage area (details on how these variables are calculated are reported in the Supporting Information). Specifically, the input concentration of all resources decreases as drainage area increases (Fig. 2d) because of the corresponding increase in water volume along the downstream direction (Fig. 2c). Note that local inputs of CPOM depend on river width and not on drainage area, therefore CPOM follows a slightly different pattern than other resources in Fig. 2d (see Supporting Information for more details). Light availability l i peaks at intermediate values of drainage area (Fig. 2e) and is lower both in upstream (due to increased shading effect of canopy in narrower reaches) and downstream reaches (due to increased water depth and subsequent limited light penetration).

We determined the equilibrium densities of each compartment of the meta-ecosystem in all reaches of the river network by finding a feasible equilibrium state (i.e., non-negative densities for all state variables (including resources) at all reaches) for system Eq. ( 1) by using a linearization method (see Supporting Information for details). In a default simulation, we investigated how the densities of functional groups change along a gradient of drainage area and compared the model predictions with the empirical patterns described in the River Continuum Concept [START_REF] Vannote | The river continuum concept[END_REF]. Here, we did not aim at testing and validating the predictions of the RCC. Rather, we aimed at building a mechanistic framework that accurately reproduces the spatial patterns described in the RCC, which we then used to explore alternative scenarios of landscape structure and resource dynamics. We also performed a sensitivity analysis (methodology and results are reported in the Supporting Information), where we varied two key hyper-parameters shaping food chains, namely predator feedback and pyramid top-heaviness (following [START_REF] Barbier | Pyramids and cascades: a synthesis of food chain functioning and stability[END_REF]).

Predator feedback expresses the ratio between consumption rates and self-regulation terms, while top-heaviness expresses the ratio between assimilation efficiency and metabolic costs.

Alternative scenarios of landscape structure and resource dynamics

We assessed the influence of assuming a complex, dendritic landscape on the spatial distribution of functional groups by running the meta-ecosystem model (Eq. ( 1)) on three alternative landscapes based on different assumptions with respect to dendritic structure and scaling of hydrological variables. In particular, we introduced a linear river channel (without branches) where hydrological variables do not scale with drainage area; a linear channel where hydrological variables scale with drainage area; and a dendritic river network where hydrological variables do not scale with drainage area. We maintained the following quantities equal in the four landscapes: number of reaches, total drainage area, water discharge at the outlet, total river length, total water volume, and total resource inputs from the terrestrial realm. Details on the construction of the alternative landscapes are provided in the Supporting Information. We then compared the spatial distributions of functional groups between the default dendritic network and the three alternative landscapes.

For the dendritic landscape, we further compared the spatial distributions and regional biomass of each functional group between the default case and two alternative scenarios of resource spatial dynamics: a "no flow" scenario, where hydrological transport of resources is neglected, and a "fast" scenario where all resources are transported at a high rate (i.e., same velocity as water). All the scenarios were based on the same set of parameters specified in Table S1.

Results

The meta-ecosystem model for dendritic river networks accurately reproduces the distributions of functional groups described in the River Continuum Concept (Fig. 3). Specifically, the density of grazers peaks at intermediate values of drainage area (Fig. 3a,f), which is mirrored by the analogous pattern for primary producers (Fig. 3g). Both of these patterns essentially follow the spatial distribution of light availability (Fig. 2e). The density of shredders is maximum in headwaters and decreases as drainage area increases (Fig. 3b,f), while filterers tend to be more homogeneously distributed across the river network (Fig. 3d,f). Conversely, the density of collectors monotonically increases with drainage area (Fig. 3c,f). Consequently, collectors and filterers are the most abundant groups in largest river reaches. In 468 reaches (12.7% of all reaches), all located at the headwaters, collectors are predicted to go extinct. Finally, the density of predators (Fig. 3e), which is proportional to the sum of their preys' densities (Fig. 3f), is relatively constant along a gradient of drainage area. The corresponding spatial patterns of resource densities predicted by the model are shown in Fig. S1. Our sensitivity analysis (see Supporting Information, Fig. S2) shows that the spatial patterns of consumers' density present the same shape with respect to the default case, while regional abundances vary substantially.

We found that all of the alternative landscapes that neglected either the dendritic structure of the river network or the scaling of hydrological variables with drainage area yielded patterns of consumer and resource densities that differ substantially with respect to the default OCN case (Fig. 4, Fig. S3). Remarkably, the linear setting without scaling of hydrological variables leads to constant patterns of consumer and resource density (a formal proof for this fact is provided in the Supporting Information). In the case of the linear channel with scaling of hydrological variables, densities of all functional groups (excluding filter feeders) tend to be much lower than in the default case. In particular, shredders' density (Fig. 4c) tends to peak for intermediate values of drainage area, and moderately decline further downstream. Collectors' density (Fig. 4d) is null in the 2% most upstream reaches, and increases downstream (with values much lower than those predicted for the default case, except towards the outlet, where values tend to be more similar).

Peak densities of grazers (Fig. 4b) and producers (Fig. 4f) are localized closer to the outlet than it was the case for the default setting. As a result, in this setting predators' density (Fig. 4a) is predicted to increase in the downstream direction. The dendritic landscape without scaling of hydrological variables leads to patterns of consumers that, in most cases, tend to decrease in the downstream direction (especially for predators, shredders, grazers and producers), while collectors' density appears highly variable in the upstream part of the network and tends to the density observed in the default case towards the outlet. In all three settings, filter feeders' density (Fig. 4e) is spatially constant and very similar to that of the default case; indeed, both the input of UPOM in the river and the consumption process by this functional group were assumed to be independent of any reach-specific variable (such as river width, which is key for CPOM input and thus shredder density, and depth, which influences producers' assimilation rate via the light limitation factor). Additional comments on these results are reported in the Supporting Information.

Concerning the effect of variations in resource flow dynamics on riverine community patterns, the absence of hydrological transport of resources has a major effect on the spatial distribution of all functional groups in the "no flow" scenario (Fig. 5). Indeed, in this scenario, the densities of all groups decrease in the downstream direction. Given that all organisms need to rely on local terrestrial inputs in this scenario, the resulting consumer patterns essentially mirror those of input resource concentration (Fig. 2d), which tend to decrease downstream due to the dilution effect (i.e., increasing water volume downstream (Fig. 2c)). This effect is particularly strong for grazers, filterers and producers (Fig. 5b,e,f), while shredders are the least impacted. Hence, functional groups feeding on fast-flowing resources (i.e., nutrients and UPOM) are more impacted than the ones feeding on slow-flowing resources (i.e., CPOM) in this scenario. For the parameter set chosen, the regional biomasses of grazers, producers and filterers are 4.8%, 6.7% and 89.1% higher in the "no flow" scenario, respectively. For these groups, the increase in density observed in upstream reaches exceeds the decrease observed downstream. Note that this result is sensitive to the type of functional response (i.e., the relationship between a consumer's intake rate and resource density [START_REF] Holling | Some characteristics of simple types of predation and parasitism[END_REF]) used in the model (see Supporting Information for additional analyses with a type II functional response). Conversely, the regional biomass of shredders is 9.4% lower, suggesting that the decrease in CPOM in the most downstream reaches has a significant effect on the regional biomass of this group. The group of collectors exhibits a totally different spatial distribution in the "no flow" scenario, with very high densities in headwaters and a decreasing density in the downstream direction (Fig. 5d). The regional biomass of collectors is more than doubled, which mirrors the amount of FPOM available in this scenario (Fig. S4). As a result of changes in prey patterns, the spatial distribution of predator density changes as well (Fig. 5a), with higher densities in upstream reaches and a regional biomass that is 54.5% higher without hydrological transport of resources.

In the "fast" scenario, increasing the rates of hydrological transport of CPOM and FPOM has major effects on the spatial distributions and overall quantity of resources available for shredders and collectors at the regional scale (Fig. 6). The amounts of CPOM and FPOM available in the whole river network are 90.8% and 99.9% lower in this scenario, respectively (Fig. S5). These changes have dramatic negative effects on the density of shredders, which is 93.6% lower at the regional scale, but shows an increasing trend in the downstream direction (Fig. 6c). We also find that shredders go extinct in 28.7% of the reaches (1058 reaches). Because FPOM is a by-product of shredders activity (Fig. 1), the sharp decrease in shredder density cascades to collectors, which go extinct in all reaches in this scenario (Fig. 6d). Consequently, the regional biomass of predators shrinks by more than half (Fig. 6a), with a spatial distribution following those of grazers (Fig. 6b) and filterers (Fig. 6e), which are unaffected by increases in the rates of hydrological transport of CPOM and FPOM.

Discussion

By performing the first formal integration of the verbal arguments of the River Continuum Concept -a milestone that has been shaping the scientific thinking in freshwater ecology over the last four decades [START_REF] Vannote | The river continuum concept[END_REF][START_REF] Doretto | The river continuum concept: lessons from the past and perspectives for the future[END_REF] -within a meta-ecosystem model, we

showed that spatial distributions and regional biomasses of major functional groups observed in stream communities are jointly shaped by the dendritic structure and scaling attributes of river networks as well as specific rates of resource flows. Neglecting any of these aspects in modelling riverine meta-ecosystems would result in very different community patterns. More generally, we showed that spatially explicit meta-ecosystem models allow understanding the interactive effects of landscape structure and resource spatial dynamics on the composition of local communities.

Model predictions did not reproduce the spatial distribution of functional groups observed in stream communities when only local ecosystem dynamics were implemented (i.e., the "no flow" scenario), which highlights the central role exerted by hydrological transport of resources in riverine ecosystems. Furthermore, the spatial distributions of functional groups were obtained without making specific assumptions on the spatial distribution of terrestrial resources (e.g., forests being more abundant at higher or intermediate elevation). For instance, in our model the spatial variation of CPOM concentration only results from the scaling of river width and water volume with drainage area in the absence of shredders, while it is modulated by shredders' assimilation and downstream transportation rates when shredders are present. Hence, we found that dendritic connectivity, hydrology-mediated resource flow and scaling of physical variables along a river network are sufficient to explain the widely observed community patterns postulated by the RCC. Moreover, a sensitivity analysis showed that the overall shape of these patterns is independent of the predator feedback and top-heaviness levels implemented in the food-web (see Supporting Information). Note that our work focuses on the spatial variations of functional group densities, not on the comparison of densities between functional groups in a given river reach, which would require a precise estimation of group-specific food-web parameters from empirical data, which is beyond the scope of our work.

We showed how the spatial distribution of functional groups crucially depends on the ratio between the feeding rate of each group and the rate at which their resources are transported across the river network. In particular, shredders are adapted to feed on CPOM, which is a slow-flowing resource. During their feeding activity, shredders produce FPOM, which is the resource on which collectors feed. We showed that, when the transport of CPOM is accelerated, which could be for instance due to the anthropogenic elimination of low-current areas or regularization of the river bed, populations of shredders are heavily depleted and this, in turn, leads to the extinction of collectors across the whole river network via a massive decrease in FPOM concentration.

Note that, in nature, part of the in-stream FPOM is likely to enter the river network directly from the terrestrial environment (through soil erosion for instance). In this case, the resulting pattern of collectors would appear as a combination of the patterns of collectors and filter feeders.

Overall, these results support the use of metacommunity and meta-ecosystem theory for ecosystem management and restoration, as suggested by Harvey et al. (2017b), by outlining the negative effect of high rates of resource transportation in river networks on several functional groups of macroinvertebrates. The meta-ecosystem approach also formalizes the relationship between the size of the particles constituting a given resource and the rate at which it is transported via hydrological flow. This relationship is based on general physical laws and could be applied to other types of resources and other freshwater organisms, such as fishes or microbes.

Our meta-ecosystem model for river networks depends on a number of technical assumptions inherited from the RCC. First, the model does not include the dispersal of living organisms, hence making the implicit assumptions that the time scale of resource transport (which is controlled by water velocity) is much shorter than that of organismal movement across reaches, and that dispersal is unbiased, with mean distances travelled that are much smaller than the extent of the river network. In particular, these assumptions allowed us to highlight the crucial effect of resource flows on the spatial patterns of riverine communities. All these hypotheses are reasonable first approximations, yet we acknowledge that more complex dispersal dynamics can have an important role in shaping riverine metacommunities [START_REF] Altermatt | Diversity in riverine metacommunities: a network perspective[END_REF][START_REF] Lowe | Is dispersal neutral?[END_REF][START_REF] Tonkin | The role of dispersal in river network metacommunities: Patterns, processes, and pathways[END_REF] and that other freshwater organisms, such as fish, can constitute important flows of energy and nutrients, both within a river network (e.g., stream-resident fishes) and across ecosystems (e.g., potamodromous or diadromous fishes).

Second, the RCC is inspired from natural, temperate-climate rivers and the predictions of the model thus apply to natural river networks that are not subject to disconnections between river reaches due to physical barriers (e.g., dams and reservoirs) or major drying events. However, a significant fraction of river networks worldwide are fragmented [START_REF] Grill | Mapping the world's free-flowing rivers[END_REF][START_REF] Belletti | More than one million barriers fragment europe's rivers[END_REF] and/or experience flow intermittence [START_REF] Allen | River ecosystem conceptual models and non-perennial rivers: A critical review[END_REF][START_REF] Messager | Global prevalence of non-perennial rivers and streams[END_REF], which is likely to disrupt the longitudinal transport of resources through the network. The effect of network fragmentation on the predictions of the RCC will depend on the location of the physical barriers or dry streams; in this respect, the use of a physically-based river landscape model as done here would enable an adequate assessment of the effects of physical barriers [START_REF] González-Ferreras | Effects of altered river network connectivity on the distribution of Salmo trutta: Insights from a metapopulation model[END_REF] and expansions/contractions of the river network [START_REF] Giezendanner | A note on the role of seasonal expansions and contractions of the flowing fluvial network on metapopulation persistence[END_REF]) on stream communities.

Recent studies outlined the need for meta-ecosystem ecology to move from a very simplified and abstract representation of ecosystems to a more realistic one (Gounand et al., 2018a;[START_REF] Guichard | Recent advances in metacommunities and meta-ecosystem theories[END_REF]. Previous theoretical studies on meta-ecosystems demonstrated that high rates of resource flow destabilize simple producer-consumer dynamics [START_REF] Marleau | Meta-ecosystem dynamics and functioning on finite spatial networks[END_REF]. However, most studies on meta-ecosystem dynamics have focused on the effect of recycling or organism movement on ecosystem stability and productivity in small spatial networks. Here, we demonstrated that meta-ecosystem dynamics and community composition are strongly influenced by the nonrandom structure of real-world landscapes. Indeed, the spatial variation of functional groups observed in river systems only emerged from the dendritic structure and scaling of physical variables characteristic of river networks, while none of the alternative landscapes neglecting either of these aspects yielded correct predictions.

While our meta-ecosystem model was designed for river networks, the general approach proposed can be adapted to other ecological systems where meta-ecosystems dynamics are important structuring processes, such as marine shorelines, coral reef systems or estuaries [START_REF] Menge | Are meta-ecosystems organized hierarchically? A model and test in rocky intertidal habitats[END_REF][START_REF] Spiecker | Reciprocal feedbacks between spatial subsidies and reserve networks in coral reef meta-ecosystems[END_REF]Gounand et al., 2018a). The application of meta-ecosystem theory to real-world landscapes crucially depends on an informed knowledge of the specific physical attributes of a landscape that constrain the direction and rates of resource flows among localities [START_REF] Polis | Food webs at the landscape level[END_REF]. Hence, the development of theoretical models, such as the one described in the present study, further outlines the need for empirical quantifications of cross-ecosystem flows of resources (e.g. [START_REF] Gounand | Cross-ecosystem carbon flows connecting ecosystems worldwide[END_REF]). Such empirical effort is central for fostering future developments of spatially explicit meta-ecosystem approaches and predicting the large-scale effects of human-induced changes in meta-ecosystem dynamics on biodiversity and ecosystem functioning. Dissolved nutrients (N) and light are the basal resources for primary producers P, while coarse (CPOM), fine (FPOM) and ultra-fine (UPOM) particulate organic matter are the basal resources for S, C and F, respectively. Predators can feed on all primary consumers, that is G, S, C and F. By feeding on CPOM, shredders produce FPOM that constitutes the main resource for collectors. The numbers within boxes correspond to the regional biomasses of the respective functional groups (total biomass density across the whole river catchment). Note that trends in panel f correspond to the patterns shown in panels a-d. Figure 4: Effect of network structure on the spatial distribution and regional biomass of functional groups in stream communities. Comparison between the spatial distribution of functional groups in a dendritic river network where river width and depth scale with drainage area (black, same as in Fig. 3) and three alternative landscapes (linear without scaling, orange; linear with scaling, violet; network without scaling, green). Solid lines represent median values calculated over reaches corresponding to each bin of Fig. 2b, while shaded areas correspond to the 2.5 th -97.5 th percentile intervals evaluated across the bins of Fig. 2b. The numbers within boxes correspond to the regional biomasses of the respective groups; percentage values within the boxes indicate the variation in total biomass relative to the default simulation. Corresponding resource patterns are shown in Fig. S3. Figure 5: Effect of hydrological transport of resources on the spatial distribution and regional biomass of functional groups in stream communities. Comparison between the spatial distributions of functional groups with (black, same as in Fig. 3) and without (purple) hydrological transport of resources in the river network. Solid lines represent median values calculated over reaches corresponding to each bin of Fig. 2b, while shaded areas correspond to the 2.5 th -97.5 th percentile intervals evaluated across the bins of Fig. 2b. The numbers within boxes correspond to the regional biomasses of the respective groups; percentage values within the boxes indicate the variation in total biomass relative to the default simulation. The corresponding resource patterns are shown in Fig. S4. Figure 6: Effect of high rates of hydrological transport of resources on the spatial distribution and regional biomass of functional groups. Comparison between the spatial distributions of functional groups for the default scenario (black, same as in Fig. 3), and the scenario with high rates of hydrological transport for all resources (blue). Solid lines represent median values calculated over reaches corresponding to each bin of Fig. 2b, while shaded areas correspond to the 2.5 th -97.5 th percentile intervals evaluated across the bins of Fig. 2b. The numbers within boxes correspond to the regional biomass of the respective groups; percentage values within the boxes indicate the variation in total biomass for the relative resource as compared with the default simulation. The corresponding resource patterns are shown in Fig. S5.

CPOM FPOM

Figure 1 :

 1 Figure1: Conceptual illustration of the meta-ecosystem model for river networks. A local stream community (grey box) is composed of a primary producer P and five functional groups of consumers (predators R, grazers G, shredders S, collectors C, and filter feeders F). The black arrows illustrate the feeding links between consumers and resources (blue box). The resources available for a given local community originate both from local terrestrial inputs (yellow ellipse on the left) and from the ecosystems situated upstream via hydrological transport (blue arrow). Dissolved nutrients (N) and light are the basal resources for primary producers P, while coarse (CPOM), fine (FPOM) and ultra-fine (UPOM) particulate organic matter are the basal resources for S, C and F, respectively. Predators can feed on all primary consumers, that is G, S, C and F. By feeding on CPOM, shredders produce FPOM that constitutes the main resource for collectors.

Figure 2

 2 Figure 2: a) The dendritic river network used in the model simulation, spanning a square of area 5625 km 2 . The river network is partitioned into 3688 reaches. Three example reaches are highlighted in black, and their corresponding drainage areas (i.e., the portion of land over which water drains into a given reach) are shown in orange; a zoom-in for the smallest of these reaches is provided at the bottom-right corner. b) Distribution of drainage area values across the 3688 reaches constituting the river network, ranging from 1 to 5625 km 2 (grouped into 16 bins). Orange solid lines display the drainage area values corresponding to the three orange areas illustrated in panel a. c-e) Variation of physical attributes of the river network variables along a gradient of drainage area: c) water volume increases with drainage area, d) concentration of terrestrial inputs decreases with drainage area, and e) light availability is maximum for intermediate values of drainage area (around 100 km 2 ). In panels c-e, solid lines represent median values calculated over reaches corresponding to each bin of panel b, while shaded areas correspond to the 2.5 th -97.5 th percentile intervals evaluated across the bins of panel b.

Figure 3 :

 3 Figure3: Spatial distributions of functional groups emerging from the meta-ecosystem model. a-d) Map representations of the distributions of the four main functional groups discussed in the RCC. e-g) Variation in density of functional group densities over drainage area. Solid lines represent median values calculated over reaches corresponding to each bin of Fig.2b, while shaded areas correspond to the 2.5 th -97.5 th percentile intervals evaluated across the bins of Fig.2b. The numbers within boxes correspond to the regional biomasses of the respective functional groups (total biomass density across the whole river catchment). Note that trends in panel f correspond to the patterns shown in panels a-d.


Lucas Baudin
email: lucas.baudin@dauphine.eu

Rida Laraki

Smooth Fictitious Play in Stochastic Games with Perturbed Payoffs and Unknown Transitions

Introduction

Fictitious play (FP) [START_REF] Brown | Iterative solution of games by fictitious play[END_REF], [START_REF] Robinson | An Iterative Method of Solving a Game[END_REF]) is the oldest and most famous independent (i.e., decentralized) learning algorithm in game theory. It is a simple behavioural strategy that agents may use to repeatedly play a fixed normal form game G. Namely, at every repetition of G, each player best replies to the empirical distribution of the uncorrelated past actions of the opponents.

When G is a finite game, it has been shown that the empirical distribution of the players' actions under FP globally converges to the set of mixed Nash equilibria, when G is zero-sum [START_REF] Brown | Iterative solution of games by fictitious play[END_REF], [START_REF] Robinson | An Iterative Method of Solving a Game[END_REF]) or identical-interest [START_REF] Monderer | Fictitious Play Property for Games with Identical Interests[END_REF]). The FP procedure was extensively studied since then with numerous additional convergence results [START_REF] Miyasawa | On the convergence of learning processes in a 2x2 non-zero-person game. Levine's Working Paper Archive[END_REF], [START_REF] Berger | Fictitious play in 2×n games[END_REF], van [START_REF] Van Strien | Fictitious play in 3×3 games: Chaos and dithering behaviour[END_REF]) and a number of generalizations [START_REF] Leslie | Generalised weakened fictitious play[END_REF], [START_REF] Benaïm | Consistency of Vanishingly Smooth Fictitious Play[END_REF] to cite a few).

Unfortunately, FP has no guarantee in terms of regret. Indeed, it may well lead to the worst possible payoff in a zero-sum game if the opponent best-replies to it at every stage. Fortunately, the smooth version of FP (call it SFP), also named stochastic or regularized fictitious play, is proved to have the no-regret property [START_REF] Fudenberg | Consistency and cautious fictitious play[END_REF], [START_REF] Benaïm | Consistency of Vanishingly Smooth Fictitious Play[END_REF]) and to converge to regularized Nash Equilibria in the classes mentioned so far [START_REF] Hofbauer | On the Global Convergence of Stochastic Fictitious Play[END_REF], [START_REF] Berger | Fictitious play in 2×n games[END_REF], [START_REF] Heliou | Learning with bandit feedback in potential games. 30[END_REF], even when G is non atomic or with continuous action sets, see 36th Conference on Neural Information Processing Systems (NeurIPS 2022). [START_REF] Hadikhanloo | Learning in nonatomic games, part I: Finite action spaces and population games[END_REF], Perrin et al. [2020]). In the online optimization and learning community, fictitious play for two player games (or correlated FP for more than two players, where a player best responds to the correlated empirical past actions of the opponents) is precisely Follow The Leader (FTL), and Smooth Fictitious Play is Follow The Regularized Leader (FTRL) (see [START_REF] Shalev-Shwartz | Online Learning and Online Convex Optimization[END_REF], [START_REF] Cesa | Potential-Based Algorithms in Online Prediction and Game Theory[END_REF], [START_REF] Sabato | Online convex optimization and no-regret learning: Algorithms, guarantees and applications[END_REF], [START_REF] Kwon | A continuous-time approach to online optimization[END_REF], Giannou et al. [2021a,b]).

How can FP and SFP be extended to learn Nash equilibria in zero-sum and identical-interest games when the stage game is not fixed but changes with time (such as the stochastic game of [START_REF] Shapley | Stochastic Games[END_REF])?

The difficulty is that, in stochastic games (SG), a state variable evolves from a period to the next one with a probability depending on the current state and actions players take. Therefore, every player, when choosing its action at any given stage, must strike a balance between its instantaneous reward (determined by the current state and the current profile of actions) and its continuation future rewards, which depend in particular on the next state. Extending the result of [START_REF] Shapley | Stochastic Games[END_REF] in the zero-sum case, [START_REF] Fink | Equilibrium in a stochastic $n$-person game[END_REF] proved that when there are finitely many players and actions (the framework of our paper), any discounted stochastic game (DSG) admits a stationary Nash equilibrium in mixed strategies (e.g. a decentralized randomized policy that depends only on the state variable).

Stationary Nash are the simplest possible equilibria that a dynamic game can have and it is our objective to construct independent learning algorithms that are provable to globally converge to them in finite zero-sum and identical-interest DSG. This is a hard problem because the payoff function of a player in a DSG is not own-payoff concave or quasi-concave when the players are restricted to their stationary strategies and thus, no gradient based method is guaranteed to converge, even to a local Nash equilibrium [START_REF] Daskalakis | The complexity of constrained min-max optimization[END_REF]). Only recently some positive results have been obtained. That is, [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF], [START_REF] Muhammed | Fictitious play in zero-sum stochastic games[END_REF], [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF] combined the Fictitious Play behavioral strategy with a Q-learning like updating rule to design a family of decentralized rules such that the sequence of independent empirical distributions of the players converge to the set of stationary Nash equilibria in the two major classes of zero-sum and identical-interest ergodic DSG.

These latest results are based on ideas of the well-known and comparatively much better understood framework of Reinforcement Learning (RL). In this setting, there is only one player (a Markov Decision Process). The advent of efficient and model-free algorithms in this context such as Qlearning [START_REF] Watkins | Learning from Delayed Rewards[END_REF]) has had a lot of impacts on RL with numerous extensions, including offline Q-learning (Kumar et al. [2020]), double Q-learning (Hasselt [2010]) and a wide range of applications [START_REF] Tai | A robot exploration strategy based on Q-learning network[END_REF], Kurin et al. [2020]). However, the convergence of Q-learning does not extend to the multiagent setting (see Wunder et al. [2010], [START_REF] Kianercy | Dynamics of Boltzmann Q learning in two-player twoaction games[END_REF]).

Unfortunately, the algorithms in [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF] and [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF] need the players to know precisely the model from the beginning (i.e., the transition and payoff functions). To avoid this drawback, our paper introduces a robust regularized version of their rules which combines Smooth FP, Q-learning-like rule (i.e., updates at every step an estimate of the continuation payoff) and empirical estimates of the unknown parameters. This leads us to a family of model-free independent learning algorithms, where the payoffs and the transitions are unknown to the players, and where the stage payoffs are imperfectly observed (e.g. randomly perturbed with a zero-mean noise). Our algorithms are proved to converge to the set of regularized stationary Nash equilibria in zero-sum and identical-interest ergodic DSG (those are a subset of stationary ε-Nash equilibria, and as ε goes to zero (corresponding to a vanishing SFP), they refine the set of stationary Nash equilibria).

Another efficient class of algorithms in RL are the Projected Gradient Methods (PGM) and their stochastic version [START_REF] Williams | Simple statistical gradient-following algorithms for connectionist reinforcement learning[END_REF], [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]). Very recently, PGM and Stochastic PGM (SPGM) have been studied by Daskalakis et al. [2020] in zero-sum stochastic games to prove that when the players play independently a PGM (or a SPGM) with different time scales,1 one can approach a best stationary equilibrium iterate2 in a finite time T (ε), whenever the stochastic game is episodic. 3 Leonardos et al. [2022] proved a similar convergence result in episodic identical-interest DSG, as soon as all the players use a PGM or SPGM. 4 Those are extremely interesting and promising convergence results. Note however that they don't cover our framework because they assume the game to have a positive probability to stop (implying that their trajectories terminate almost surely in finite time) while in our ergodic setting, all the trajectories are infinite and the game never terminates. So, our results and those of Daskalakis et al. [2020], [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF] are non-comparable and complementary: they prove the existence of a best iterate approximation, we prove a time-average convergence, we use different algorithms and tools, and orthogonal assumptions (episodic vs ergodic).

The counterpart of FP in continuous time is best-response dynamics [START_REF] Matsui | Best response dynamics and socially stable strategies[END_REF], [START_REF] Harris | On the Rate of Convergence of Continuous-Time Fictitious Play[END_REF]) and modern proofs for the convergence of SFP rely on the convergence of the continuous model combined with stochastic approximation techniques (see for instance [START_REF] Benaïm | Stochastic Approximations and Differential Inclusions[END_REF], [START_REF] Benaïm | Consistency of Vanishingly Smooth Fictitious Play[END_REF], [START_REF] Hadikhanloo | Learning in nonatomic games, part I: Finite action spaces and population games[END_REF]). We follow a similar approach: we prove the convergence to regularized stationary Nash equilibria of an associated smooth continuous-time dynamics and, using some advanced stochastic approximation tools, deduce the convergence of our discrete-time rules. This is an important technical difference with Daskalakis et al. [2020], [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF] who can prove their results directly in discrete time, with an explicit finite convergence bound T (ε). Our stochastic approximation methods do not provide us with a convergence rate, which is an open problem for SFP and FP even in the classical setting.

Contributions

• We introduce a family of independent learning algorithms in stochastic games that are model-free (unknown transitions and imperfect observation of the stage payoffs); • We identify the corresponding smooth continuous-time dynamics, and show it globally converges to regularized stationary Nash equilibria in identical-interest and zero-sum DSG; • From the continuous-time convergence, we deduce that in our two classes of DSG, the uncorrelated empirical frequencies of actions generated by our discrete-time algorithms almost surely globally converge to the set of regularized stationary Nash equilibria.

Outline Section 2 gives the main definitions and notations of the paper. Section 3 introduces the smooth fictitious play procedure together with the continuation payoff updating rule. In order to prove the convergence of the discrete time algorithm, and also for its own sake, a smooth best-response dynamics with a simple continuation payoff up-dating rule dynamics is described in Section 4 followed with sketches of proofs of both our discrete-time and continuous-time systems. Section 5 describes in more details the related work. The appendix contains the detailed proofs. Appendix C explains an example and presents empirical results of our algorithm.

Preliminaries

Strategic situations where several agents interact, get rewards and modify an environment can be modeled as Stochastic Games (SG). In our settings, the number of players, actions and possible states are finite. We consider players that are interested in the so-called discounted reward on an infinite horizon, that is players strike a balance between instantaneous rewards and future ones.

Stochastic games SG are tuples G = (S, I, A, {r i s } i∈I,s∈S , {P s } s∈S )
where S is the state space (a finite set), I is the finite set of players, A i is the finite action set of player i, A := Π i∈I A i is the set of action profiles, r i s (•) : A → R is the stage reward of player i, and P s (•) : A → ∆ (S) is the transition probability map (where ∆ (S) is the set of probability distributions on S).

Main Restrictions

We are interested in two classes of games: zero-sum stochastic games are two-players SG where r 1 s = -r 2 s for every s and team stochastic games are such that the payoff functions of the players differ only by a constant (there is r s (•) : A → R such that for every i and s, r i s (•) = r s (•) + c i for a constant c i ). A special case is identical-interest SG where all payoff functions are equal (r i s (•) = r s (•) for all i). An SG is ergodic if there exists T ∈ N such that for any sequence of actions of length T , the probability to reach any state s ′ starting from any state s is positive.

Game Form A stochastic game is played as follows: starting from an initial state s 0 , at every step n ∈ N, every player i choose an action a i n given the history of play and the current state s n . The next state s n+1 is drawn from distribution P sn (a n ).

Discounted Payoffs We suppose that every player i is interested in maximizing its discounted payoff, that is the expectancy of k∈N δ k r i s k (a k ) where δ ∈ (0, 1) is the discount factor.

Strategies A behavioral strategy σ i for player i is a mapping associating with each stage n ∈ N, history h n ∈ (S × A) n and current state s, a mixed action

x i n = σ i (n, h n , s) in ∆(A i ).
The behavioral strategy is pure if its image is always in A i . A stationary strategy of player i is the simplest of behavioral strategies. It depends only on the current state s but not on the period n nor on the past history h n . As such, a stationary strategy can be identified with an element of ∆(A i ) S (a mixed action per state interpreted as: whenever the state is s, i plays randomly according to x i s ). The set of stationary strategy profiles is Π i∈I ∆(A i ) S . Set X = Π i∈I ∆ A i , so a stationary profile is an element of X S . For y i ∈ (∆(A i )) S and x ∈ Π i∈I ∆(A i ) S , we denote by (y i , x -i ) the stationary profile where i changes its strategy from x i to y i . A stationary profile x ∈ X S is a Nash equilibrium if and only if no player has a profitable behavioral deviation. [START_REF] Fink | Equilibrium in a stochastic $n$-person game[END_REF] proved the existence of stationary Nash equilibria in every finite DSG and that it is sufficient to check pure stationary deviations.

Regularizer In this paper, we are interested in exploratory algorithms, which may classically be generated by the use of some steep concave regularizer. This regularizer is added to the payoff functions r i s and can be given several interpretations (see [START_REF] Fudenberg | The Theory of Learning in Games[END_REF], [START_REF] Hofbauer | On the Global Convergence of Stochastic Fictitious Play[END_REF] for details): it models the uncertainty of the payoff caused by the "trembling hand" of players or is a way to generate strict incentives to explore all the actions. Formally player i maximizes a perturbation of its payoff function r i s + ϵh i under the following hypotheses: h i : X → R + , strictly concave in x i , C 1 on the interior,

lim x i →∂∆(A i ) ∥∇ x i h i (x) ∥ = +∞ and ϵ > 0 (H1)
In this paper we study the convergence of some discrete and continuous time systems to regularized Stationary Nash equilibria, parameterized by the regularizers (h i ) i∈I and parameter ϵ. Definition 2.1. Regularized Stationary Nash Equilibria of a DSG are the stationary Nash equilibria of the DSG with the perturbed payoff functions r i s + ϵh i , i ∈ I.

Similarly to stationary Nash equilibria, there exists at least one regularized stationary Nash equilibria [START_REF] Takahashi | Equilibrium points of stochastic non-cooperative n -person games[END_REF].

Remark. In identical-interest (resp. zero-sum) SG, we suppose all players take the same (resp. the opposite) regularizer h so as the regularized payoff functions remain identical. In this context, one can take for example a separable regularizer function which is equal to the sum (resp. the difference) of concave functions depending only on x i . Remark. If in a profile, all actions are β-optimal with respect to discounted payoff based on functions r i s + ϵh i (i.e., deviations are at most β profitable), then this profile is called a β-regularized equilibria.

3 Smooth Fictitious Play in Stochastic Games (Discrete-Time Algorithms)

Fictitious Play in Repeated Games Introduced by [START_REF] Robinson | An Iterative Method of Solving a Game[END_REF] and [START_REF] Brown | Iterative solution of games by fictitious play[END_REF], fictitious play is a decentralized behavioral strategy to repeatedly play a fixed normal form game. At every step n, every player i chooses an action a i n+1 that is a best response to the past empirical average action of other players:

a i n+1 must maximize r i •, x -i n where x -i n = (x j n ) j̸ =i and ∀j ∈ I, x j n = 1 n+1 n k=0 a j k .
A famous variation of this procedure is smooth fictitious play [START_REF] Fudenberg | Consistency and cautious fictitious play[END_REF]) where players choose their action according to a regularized payoff function. Formally, a player i draws an action according to a distribution that maximizes r i •,

x -i n + ϵh i •, x -i n
(with h i and ϵ defined in the previous section). An interesting property of such a procedure is that with suitable property on ϵ, it has no regret (up to ϵ), meaning that if played unilaterally by a player i, other players can not trick i into using a suboptimal action. This is due to the randomness of the action choice: the distribution assigns positive probability to every action, so a player's behavior remains unpredictable.

In this section, we extend the smooth FP procedure to stochastic games. This builds upon the recent extension of FP in [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF], [START_REF] Muhammed | Fictitious play in zero-sum stochastic games[END_REF], [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF] (see Section 5 for more details). Indeed, it is not easy to derive the convergence to a regularized equilibrium by applying directly the definition of smooth fictitious play to the discounted stochastic game in which the players are restricted to play in stationary strategies, because the payoff function in this game is non-linear (nor is it concave, or quasi-concave) with respect to a player stationary strategy. To overcome this difficulty, and following the idea in [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF], we update two sets of variables for every state: one concerns the uncorrelated empirical actions and the other is an estimate of the continuation payoffs (i.e., payoffs that players can anticipate to achieve if that state is reached). These continuation payoffs are used as a parameter in an auxiliary game, often called the Shapley operator.

Auxiliary game Following numerous authors including [START_REF] Shapley | Stochastic Games[END_REF], we define a so-called auxiliary game parameterized by a family of vector u i ∈ R S (one for every player i) and a state s. It is a one-shot game every player i has the same action set as in G but gets a one-time payoff of:

f i s,u i (a s ) := (1 -δ)r i s (a s ) + δ s ′ ∈S P ss ′ (a s )u i s ′ (1)
Function f i s,u i is extended to mixed action profiles.

Smooth Best-Response While playing the stochastic games, players maintain a set of continuation payoffs u i that are used to choose their actions. Given that the current state is s, for a player i, its action is drawn from the distribution that is the smooth best-response in the auxiliary game with respect to empirical action profile x s , that is:

sbr i s,u i (x s ) := arg max y i ∈∆(A i ) f i s,u i y i , x -i s + ϵh i y i , x -i s (2)
This is well and uniquely-defined because of the strict concavity of h i .

Example. If the regularizer h i is taken to be the Shannon entropy, that is:

h i (x s ) = - j∈I a j ∈A j
x j s (a j ) log x j s (a j )

Then the smooth best-response function is the logit function:

sbr i s,u i (x s ) (a i ) := exp(ϵ -1 f i s,u i a i , x -i s ) b i ∈A i exp(ϵ -1 f i s,u i b i , x -i s )
See [START_REF] Mertikopoulos | Learning in Games via Reinforcement and Regularization[END_REF] for other examples of regularizers.

Smooth fictitious play with known transitions and deterministic payoff To extend smooth fictitious play to stochastic games, we use two sets of variables. The x s,n variable is the distribution of empirical actions for every state s prior to step n. The other variable u i s,n+1 can be interpreted as the continuation payoffs at s and is defined as the time average of the regularized payoffs up to step n. This leads to the following system:

                 u i s,n+1 = 1 n + 1 n k=0 f i s,u i k x s,k + ϵh i x s,k x s,n = 1 s ♯ n n k=0 1 s=s k a k a i n+1 ∼ sbr i s,u i n x s,n
where s ♯ n is the number of times s was reached, that is n k=0 1 s=s k and every action a k is embedded into the Euclidean space containing ∆ A i . The interpretation is simple: at each stage t, each player best replies to the belief that the other players will play according to the past uncorrelated empirical distribution of actions, and that the future continuation payoffs are equal to the time average of the past estimated perturbed payoffs, calculated using the past empirical frequencies of actions.

This system can be generalized and rewritten in an incremental fashion. It is our first main system:

             u i s,n+1 -u i s,n = β n + 1 f i s,u i n x s,n + ϵh i x s,n -u i s,n x s,n+1 -x s,n = 1 s=sn s ♯ n a i n -x s,n a i n+1 ∼ sbr i s,u i n x s,n (SFP)
where β = 1 corresponds to the equation above.

Remark. If there is only one state (e.g the classical repeated game sitting), this procedure is exactly standard smooth fictitious play since the smooth best-response does not depend on the value of u i s,n .

Smooth fictitious play with unknown transitions and perturbed payoff We suppose now that payoff functions are not known, and that each stage payoff is observed with some zero-mean noise which follows a distribution that may depend on the history, the current state and actions taken by the other players. Therefore, at step n, player i gets a random reward R n i that is drawn according to a distribution determined by actions a n and current state s n whose expectancy is r i sn (a n ) and bounded variance conditionally on the history (see Appendix B for details). We also suppose that transitions are not known. Therefore, both transitions and expected payoff may be empirically estimated as follows:

         Pss ′ ,n (a) = n k=0 1 s k =s∧a k =a 1 s k+1 =s ′ n k=0 1 s k =s∧a k =a ri s,n (a) = n k=0 1 s k =s∧a k =a R k i n k=0 1 s k =s∧a k =a (MFP.1)
Consequently, we define the estimated auxiliary payoff using these two estimators:

f i s,u i (x s ) := (1 -δ)r i s,n (x s ) + δ s ′ ∈S Pss ′ ,n (x s ) u i s ′ (MFP.2)
Now we define a model-free version of smooth fictitious play, similar to SFP but using estimators:

             u i s,n+1 -u i s,n = β n + 1 f i s,u i n x s,n + ϵh i x s,n -u i s,n x s,n+1 -x s,n = 1 s=sn s ♯ n a i n -x s,n a i n+1 ∼ sbr i s,u i n x s,n (MFP)
where sbr i s,• is defined relatively to f i s,• . Theorem 3.1 (Convergence in identical-interest ergodic DSG). In a identical-interest ergodic discounted stochastic game, if all players follow SFP (resp. MFP), their empirical actions x s,n converge almost surely to the set of regularized stationary Nash equilibria and their expected vector of continuation payoffs u i s,n converges to the optimal continuation payoff of limiting equilibrium set.

This theorem means that even if the trajectories of SFP (resp. MFP) cycle between several stationary equilibria, they all share the same optimal continuation payoff vector. Theorem 3.2 (Convergence in zero-sum ergodic DSG). In a zero-sum ergodic discounted stochastic game, if all the players follow SFP (resp. MFP) with the same initial values, their empirical actions x s,n converge almost surely to the set of Dβ-regularized Nash equilibria (where D > 0 is a constant that only depends on G) and their expected vector of continuation payoffs u i s,n converges to the corresponding continuation payoff.

To prove these results, we are going to define in the next section an associated smooth best-response continuous time counterpart to our independent learning algorithms. The proof of Theorems 3.1 and 3.2 are sketched at the end of Section 4 and fully detailed in Appendix B. Remark 1. Theorem 3.2 suggests that it is possible to use a doubling-trick mechanism to converge to the set of 0-regularized stationary Nash equilibria. Indeed, players can compute the duality gap (see Appendix A for a definition) and decide to reduce the update rate β every time it is below a certain threshold. This is a standard trick to achieve no-regret in reinforcement learning. Remark 2. Our rules suppose that each player observes the other player's past actions. This is an important difference with the algorithms developed in Daskalakis et al. [2020], [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF] where a player observes only its own actions. On the other hand, they do not prove a time-average convergence of their trajectories, while we do, but they prove a best iterate convergence that we do not prove.

Smooth Best-Response in Stochastic Games (Continuous-Time Dynamics)

Continuous counterpart of discrete time systems The continuous-time counterpart of fictitious play is best-response dynamics of [START_REF] Matsui | Best response dynamics and socially stable strategies[END_REF], [START_REF] Harris | On the Rate of Convergence of Continuous-Time Fictitious Play[END_REF]. Together with the theory of stochastic approximations [START_REF] Benaïm | Stochastic Approximations and Differential Inclusions[END_REF]), proofs of convergence of best-response are key ingredients to modern proofs of convergence of fictitious play and smooth FP in repeated games [START_REF] Benaïm | Stochastic Approximations and Differential Inclusions[END_REF], [START_REF] Benaïm | Consistency of Vanishingly Smooth Fictitious Play[END_REF]) and more recently in proving the convergence of FP in stochastic games [START_REF] Muhammed | Fictitious play in zero-sum stochastic games[END_REF], [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF]).

In this section, we define the continuous-time counterpart to our discrete-time algorithms. This leads us to a regularized version of the continuous-time best-response dynamics in [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF], [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF]. Similarly to the previous section, our continuous time system has two sets of variables u i s (t) which represents the belief about the continuation payoffs and x s (t) which corresponds to belief about the players actions. Our main continuous-time system is:

         ui s = β(t) f i s,u i (t) (x s (t)) + ϵh i (x s (t)) -u i s (t) ẋs = α s (t) sbr i s,u i (t) (x s (t)) -x s (t) α s (t) ∈ [α -, 1] (SBRD)
Remark 3. In zero-sum or identical-interest games, if players have the same initial conditions, then we can omit the superscript i in u i s as they are equal.

Update rates Profile x s (t) evolve towards the smooth best-response at a rate of α s (t). Variable α s (t) corresponds to the frequency at which a state s is visited. The fact that it is bounded below by α -> 0 is a mathematically convenient way to exploit the ergodicity of the stochastic game. Remark 4. Note that this is different from the model of continuous-time stochastic game outlined for instance by [START_REF] Neyman | Continuous-time stochastic games[END_REF]. In this paper, we are interested in using the theory of stochastic approximation, therefore discrete-time and continuous-time system are related through an exponential change of variable in time. As a consequence, we can consider that the state occupation is averaged, which is not the case in some other studies in stochastic games.

Regarding continuation payoffs, they evolve towards perturbed payoffs at a rate of β(t), which is typically oblivious, that is only determined by time t and state s independently of the value of other variables. We make the following assumptions that guarantee that the update rates are not too small:

β(t) ≥ 0 and β is decreasing ∞ 0 β(u)du = +∞ (H2)
Remark 5. With β(t) = 1 t+1 , the system is close to the one outlined in [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF]. In this case, continuation payoffs are updated with a slower timescale. Indeed, having two timescales may be useful to establish convergence or apply stochastic approximation theorems. In our paper, we do not make such an assumption, allowing the possibility that the two variables are on a single timescale.

Model-Free System

In order to also study the convergence of the model-free version of our procedure, we also define the smooth best-response dynamics when the model is progressively learned. The estimators are defined as follows for every pure profile b ∈ A:

       Ṗss ′ (b) = α s (t)a (t) (b) P ss ′ (b) -Pss ′ (b) (t) ṙi s (b) = α s (t)a (t) (b) r i s (b) -ri s (b) (t) a i (t) = sbr i s,u i (t) (x s (t)) (3)
where a (t) := Π i∈I a i (t) is the profile of selected actions and a (t) (b) is the joint probability to select pure profile b at time t.

Then, f i s,• is defined as in (MFP.2) and we obtain a system MBRD by using estimators in SBRD.

Solutions Systems SBRD and MBRD are differential inclusions. A general theory can be found in [START_REF] Aubin | Differential Inclusions: Set-Valued Maps and Viability Theory[END_REF]. Here, the right-hand side of each system forms a closed, set-valued map (because there are several possible values for α s (t)) whose images are convex and uniformly bounded. Under these assumptions, it is known that such systems admit (typically non unique) solutions. Theorem 4.1. Under hypothesis H2, MBRD converges to the set of regularized Nash equilibria in identical-interest stochastic games. Furthermore, if β ⋆ ≥ lim sup t→∞ β(t), then MBRD converges to the set of Dβ ⋆ -regularized stationary Nash equilibria in zero-sum games, where D is a positive constant that depends only on the game G. In particular if β(t) goes to 0 then MBRD converges to the set of regularized stationary Nash equilibria.

The convergence of MBRD is helpful to characterize the limit set of discrete-time smooth fictitious play systems. Indeed, they are contained in the internally chain transitive sets (see Appendix B for a definition) of MBRD using the theory of stochastic approximations. Below, we sketch the proof of continuous-time Theorem 4.1. The discrete-time counterpart Theorem 3.1 and Theorem 3.2 can be deduced in a similar fashion as [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF], although there are several technical subtleties. Complete proofs of these results are in Appendix A.

Sketch of the proof of Theorem 4.1. The proofs for both class of games proceed quite differently: this is not surprising since there are no (to the best of our knowledge) unified convergence proof of even simple FP in potential and zero-sum games.

For identical-interest stochastic games, the key point is showing that the gap between (estimated) auxiliary payoffs f i s,u i (t) (x s (t)) and the auxiliary values u i s (t) is narrowing. Technically, the difference is bounded below by an (absolutely) integrable function. Since it is the differential of u i s (t) and that u i s (t) it, which implies that u i s (t) converges and then that f i s,u i (t) (x s (t)) converges and their limits are necessarily equal. Then, a study of the behavior of actions shows that they belong to the set of regularized equilibria, otherwise f i s,u i (t) (x s (t)) could not converge (based on Lipschitz properties of all these quantities).

Regarding zero-sum stochastic games, the first part of the proof studies a rather standard quantity called the duality gap. It goes to 0 (or at least near 0), which implies first that the value of the auxiliary game is mostly reached by the auxiliary payoffs f i s,u i (t) (x s (t)) and second, that the minmax strategies of the auxiliary game is learned by the players. Then, comparing relative speed of auxiliary values in all states leads to the convergence of values u i s and of other variables.

5 More on the Related Works Single-Player Markov Decision Process [START_REF] Watkins | Learning from Delayed Rewards[END_REF] introduced Q-learning to solve Markov Decision Process with guarantees when the environment is stationary. The central idea is to estimate the so-called Q-values for every state-action pairs that represents the continuation payoff using the following scheme:

Q(s, a) ← Q(s, a) + β(r s (a) + δ max a ′ Q(s ′ , a ′ ) -Q(s, a)) ( 4
)
where s is the current state, a the action played, s ′ the next state. Q-values for other state-action pairs are left unchanged. This simple rule does not need any information about the transitions.

Since then, numerous extensions have been introduced, for instance double Q-learning (Hasselt [2010]) or Q-learning with deep learning [START_REF] Hado Van Hasselt | Deep reinforcement learning with double qlearning[END_REF]). However, in multiagent systems, from the one player's point of view the environment is not stationary because it comprises other players that are also using some non stationary learning rule. Therefore, there is no guarantee and indeed it may not converge to the set of Nash equilibria (Wunder et al. [2010], [START_REF] Calvano | Artificial Intelligence, Algorithmic Pricing, and Collusion[END_REF].

Our paper is a way to overcome this difficulty with a different updating rule for estimating the continuation payoffs. This updating rule is closer to Expected Sarsa [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]) where the continuation values are moved towards the expectation of future payoffs.

Nash Q-learning [START_REF] Hu | Nash q-learning for general-sum stochastic games[END_REF] is another way to learn Nash equilibria based on Qlearning, supposing that players have access to a Nash equilibrium of the auxiliary game parameterized by the Q-values. In our paper, we do not suppose that the players have access to such an oracle.

Fictitious Play Algorithms in Stochastic Games [START_REF] Brown | Iterative solution of games by fictitious play[END_REF] and [START_REF] Robinson | An Iterative Method of Solving a Game[END_REF] introduced fictitious play to learn the value of repeated fixed zero-sum game in a turn based fashion (see also [START_REF] Berger | Brown's original fictitious play[END_REF]). It was later studied in the usual simultaneous updating by [START_REF] Fudenberg | The Theory of Learning in Games[END_REF] with numerous extensions, including smooth fictitious play. However, none of these procedures converge to equilibria in the repetition of the zero-sum normal form game Γ induced by the stationary strategies of our DSG, mainly because the payoff functions of Γ, on a player own-strategy, are neither linear, nor concave [START_REF] Daskalakis | The complexity of constrained min-max optimization[END_REF] and so new idea are necessary to overcome this difficulty. [START_REF] Vrieze | Fictitious play applied to sequences of games and discounted stochastic games[END_REF] studied fictitious play in the context of varying stage games. However, similarly to the original fictitious play of Robinson and Robinson, it is designed as a way to compute the value of a zero-sum game, and it is not a behavioral strategy (i.e., there is no current state, the goal is to compute the minmax of a series of matrices). More recently, [START_REF] Leslie | Best-response dynamics in zero-sum stochastic games[END_REF] proposed a best-response dynamics in continuous time which converges to stationary Nash equilibria in zero-sum DSG. The convergence was extended to discrete-time procedures and to identical-interest discounted stochastic games by [START_REF] Baudin | Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games[END_REF] Smooth or Regularized Learning Learning using regularizers is a widespread technique in machine learning, for instance with "follow the perturbed leader" [START_REF] Cesa | Prediction, Learning, and Games[END_REF]) or so-called stochastic fictitious play [START_REF] Fudenberg | Consistency and cautious fictitious play[END_REF]). It makes it possible to have no-regret [START_REF] Perchet | Approachability, regret and calibration: Implications and equivalences[END_REF][START_REF] Shalev-Shwartz | Online Learning and Online Convex Optimization[END_REF]) and from another point of view it is a substitute to the simpler ϵ-greedy exploration scheme in reinforcement learning [START_REF] Sutton | Reinforcement Learning: An Introduction[END_REF]).

Optimistic Gradient Descent/Ascent [START_REF] Wei | Last-iterate Convergence of Decentralized Optimistic Gradient Descent/Ascent in Infinite-horizon Competitive Markov Games[END_REF] recently proposed another algorithm to learn in stochastic games. Similarly to ours, their method is rational in the sense that, if the opponent is playing some arbitrarily stationary strategy, then the agent learns to react optimally. However, this is restricted to zero-sum games and is a procedure quite different from ours, it uses projected gradient ascent/descent instead of fictitious-play like action selection.

Local vs Global Our algorithms lead to all the players globally optimizing. Other algorithms such as COMA [START_REF] Foerster | Counterfactual multi-agent policy gradients[END_REF]) lead to local optima, and so belong to another line of work.

Conclusion

We defined a number of decentralized continuous and discrete time systems with exploration to learn its own payoff function, the transition probability function, and which converge to regularized (and so ε-approximate) stationary Nash equilibria in discounted stochastic games. This is an area that has not been widely studied and as such, a number of questions remain to be answered-either regarding our systems or regarding other recent decentralized learning algorithms in stochastic games.

First, the theory of stochastic approximations gives no clue about the speed of convergence of the corresponding discrete-time algorithms. Moreover, even in continuous time, this question remains open even for potential games [START_REF] Harris | On the Rate of Convergence of Continuous-Time Fictitious Play[END_REF]).Therefore, an interesting question in both the stochastic and the repeated game setting is, which guarantee FP and SFP have regarding the rate of convergence? Second, and inspired by the convergence of vanishing fictitious play to equilibria in classical repeated games (see [START_REF] Benaïm | Consistency of Vanishingly Smooth Fictitious Play[END_REF], [START_REF] Hadikhanloo | Learning in nonatomic games, part I: Finite action spaces and population games[END_REF]), it would be of interest to study a vanishing version of our algorithm with a parameter ε(t) that goes to zero in a way that guarantees the convergence to (exact) stationary Nash equilibria of the DSG.

Third, convergence of FP or SFP in repeated game extends to infinite action games and to non-atomic games. It would be of interest to relax our finiteness assumption, at least for the action sets.

A challenging open problem is the design of independent learning algorithms that converge to stationary equilibria in ergodic zero-sum and identical-interest DSG without the knowledge of the other player's past actions. The Projected Gradient Methods of Daskalakis et al. [2020], [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF] have this minimal information property, but their results are proved only episodic SG.

Last but not least, obtaining a last-iterate convergence instead of a time average convergence (as in our paper) or best iterate convergence (as in Daskalakis et al. [2020], [START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF]) is another interesting direction. However, it is known that there is no necessary last-iterate convergence with smooth fictitious play in the non-stochastic case (Giannou et al. [2021a]). Therefore, we can imagine that this would require a different algorithm but we emphasize that algorithms without last-iterate convergence such as ours are still of interest because of their behavioral, micro-economics and experimental foundations [START_REF] Fudenberg | The Theory of Learning in Games[END_REF], [START_REF] Hofbauer | On the Global Convergence of Stochastic Fictitious Play[END_REF]).

  . These discrete-time and continuous-time systems are close to our systems, but they only work in model-based settings with no perturbation, whereas ours is suited for model-free settings and uses smooth best-responses.[START_REF] Muhammed | Fictitious play in zero-sum stochastic games[END_REF] introduced a fictitious play algorithm for zero-sum stochastic games in discrete time with a continuation updating rule close to Q-learning in a model-based or model-free setting: both players update a Q-table with an entry per state-action pair whereas our continuation payoff are only indexed by the state. Compared to our work, algorithm of Sayin et al.[2020] does not estimate the model explicitly, which may or may not be an advantage depending on the context. Furthermore, their action selection is not smooth-the selected action must maximize the Q-values with an ϵ-greedy strategy in the model-free case. In contrast, our work uses smooth best-responses, whose main advantage is that in the non-stochastic setting, this is known to have no-regret properties, while any ϵ-deterministic procedures are known to have regret. Future work could formalize this notion of regret in the context of stochastic games.Finally, Sayin et al. [2020] algorithm and proof have been established in the zero-sum case while we consider both zero-sum and identical-interest stochastic games.[START_REF] Muhammed | Fictitious play in markov games with single controller[END_REF] very recently introduced another Fictitious Play with provable convergence in identical-interest stochastic games with single controller, which is a restriction we do not make in this paper.

Daskalakis et al. [2020] proved that their convergence result fails if the players have the same time scale.

I.e., there is t ∈ {1, ..., T (ε)} such that (in expectation in case of Stochastic PGM), the players play an ϵ stationary Nash equilibrium at time t.

A SG is episodic if, at every stage, there is a positive probability that the game stops.

[START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF] proved their results in the larger class of potential stochastic games (they called Markov Potential Games). In such a class of games, the players are divided into two categories: either they do not influence the transition or they have the same payoff function up to a constant (see Holler[2020],[START_REF] Leonardos | Global convergence of multi-agent policy gradient in markov potential games[END_REF]). This last class is slightly larger than identical-interest SG is called Team Stochastic Games in Holler[2020]. All our convergence results extend to team stochastic games. The class of potential SG where the transitions are independent of the player's actions is easy to solve, just let each player uses a smooth FP myopically state per state independently.

Michael Wunder, Michael Littman, and Monica Babes. Classes of multiagent q-learning dynamics with epsilon-greedy exploration. In Proceedings of the 27th International Conference on International Conference on Machine Learning, ICML'10, pages 1167-1174, Madison, WI, USA, June 2010. Omnipress. ISBN 978-1-60558-907-7.

Checklist

1. For all authors...