Pierre Couteron
email: pierre.couteron@ird.fr.

Conservative or dissipative? Two distinct processes for spatial pattern emergence

In PNAS, Siteur et al. [START_REF] Siteur | Phase-separation physics underlies new theory for the resilience of patchy ecosystems[END_REF] recently contributed to the topic of spatial patterns in ecosystems. This is of relevance since ecologists have long recognized that ecological variables (e.g., plant biomass, species abundances, etc.) display uneven distributions in space-i.e., patterns-that need to be accounted for and interpreted [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF][START_REF] Rietkerk | Regular pattern formation in real ecosystems[END_REF]. Although ecologists have historically favored external factors, such as preexisting substrate heterogeneity or disturbance as potential explanations, we have witnessed an increased awareness of the potential for self-organizing processes to account for various pattern types, some of them amazing (Fig. 1). Self-organized spatial patterns have now been evidenced or hypothesized in a variety of ecosystems. Yet vegetation in drylands has become paradigmatic for the introduction of pattern formation approach (sensu 4, 5) in the study of ecosystems. Notably, the most impressive patterns displaying both spatial periodicity and stability in time [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF][START_REF] Rietkerk | Regular pattern formation in real ecosystems[END_REF] are observable worldwide across drylands in relation to predictable ranges of environmental conditions [START_REF] Deblauwe | The global biogeography of semi-arid periodic vegetation patterns[END_REF]. Spatially periodic patterns inspired a first contribution that introduced in vegetation ecology methods from the field of morphogenesis in physical and chemical systems [START_REF] Lefever | On the origin of tiger bush[END_REF]. This seminal approach proposed an integrodifferential modeling of plant biomass dynamics using spatial kernels to embody plant-plant interactions, whether positive (enhancement of scarce water resource by shading) or negative (resource preemption by roots). Pattern formation through spatial-symmetry breaking instability directly relates to relative strengths and spatial scales of these two types of interactions, a necessary condition being that positive and negative feedbacks predominate at short and longer ranges, respectively [further referred to as scale-dependent-feedbacks (SDF)] (1). In this approach, and deriving subsequent works (e.g., refs. 8-10), the fast dynamics of the limiting resource, namely water, remained implicit. In parallel, the subject inspired theoretical efforts that revisited and adapted the classical approach of morphogenesis from A. Turing [START_REF] Rietkerk | Regular pattern formation in real ecosystems[END_REF][START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF] and addressed the interacting dynamics of vegetation biomass and water (explicitly considered) through systems of reaction-diffusion equations (e.g., refs. 12 and 13). A comprehensive model integrating most relevant aspects of intertwined water-vegetation dynamics in drylands is detailed in refs. 4 and 5, while several specific contributions relate to it (e.g., refs. 14 and 15). Quickly too, patterns outside the arid lands alone were approached in ecosystems as varied as peat bogs, mussel beds, etc. [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF][START_REF] Rietkerk | Regular pattern formation in real ecosystems[END_REF].

Despite that archetypal, strikingly regular spatial patterns inspired most of the aforementioned efforts, it is obvious that a dominant share of observable real-world patterns are irregular. Understanding patterns, whether regular or irregular, not only matters for grasping their role in ecosystem functioning but is also hypothesized as a way to get indicators, or "early warning signals" of possible ecosystem collapse under external stressors, such as drought, overgrazing, etc. [START_REF] Rietkerk | Self-organized patchiness and catastrophic shifts in ecosystems[END_REF][START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Rietkerk | Evasion of tipping in complex systems through spatial pattern formation[END_REF]. This question remains timely although several recent contributions have rather underlined how the relations between patterns properties and the dynamics in progress are often complex [START_REF] Sherratt | Using wavelength and slope to infer the historical origin of semiarid vegetation bands[END_REF][START_REF] Zelnik | Gradual regime shifts in fairy circles[END_REF][START_REF] Rietkerk | Evasion of tipping in complex systems through spatial pattern formation[END_REF][START_REF] Sankaran | Clustering and correlations: Inferring resilience from spatial patterns in ecosystems[END_REF]. Irregular patterns are of different sorts and have been interpreted following different theoretical lines. "Localized structures" have been frequently highlighted as additional outcomes of models aimed at explaining regular patterns [START_REF] Couteron | Plant clonal morphologies and spatial patterns as self-organized responses to resource-limited environments[END_REF][START_REF] Bordeu | Self-replication of localized vegetation patches in scarce environments[END_REF][START_REF] Zelnik | Gradual regime shifts in fairy circles[END_REF]. In spite of irregular spatial distribution, localized structures display as signature a dominant patch size. Scale-free spatial patterns (14-16) have been studied and hypothesized as dynamical outcomes in ecological systems displaying strong positive feedbacks and devoid of negative feedbacks. This precludes regular pattern formation while enabling patchiness through long-lasting phase coexistence [START_REF] Siteur | Phase-separation physics underlies new theory for the resilience of patchy ecosystems[END_REF]. In these systems, patch-size distribution is expected to follow power-laws [START_REF] Hardenberg | Periodic versus scale-free patterns in dryland vegetation[END_REF][START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF][START_REF] Sankaran | Clustering and correlations: Inferring resilience from spatial patterns in ecosystems[END_REF]. Considering the variety, the complexity, and the importance of the subject, it is obvious that the questions linked to nature, resilience, and stability of irregular patterns deserves increased attention.

In PNAS, Siteur et al. ( 1) highlighted a particular subset of patterning processes, which may result in irregular patterns having specific properties in terms of patch-size distribution and patch dynamics. Such process has been so far largely overlooked in ecology. It characteristically features the interplay between two forms of a given substance. One form is mobile and the other sessile, while the two corresponding quantities sum to a fixed amount assuming neither creation nor leakage in the whole system. This contrasts with most published models addressing patterning processes that are referred to as "dissipative" for featuring a critical substance flowing (explicitly or not) through the system. For instance, existing dryland models entail (or invoke) water dynamics that vegetation either consume or concentrate/preserve from evaporation or runoff.

From this premise, the authors propose a generic framework that may accommodate patterns resulting from phenomena as diverse as plant-nutrient cycling, picky grazing by herbivores, or sand movement in vegetated dune fields, with likely applicability to additional real-world examples. As idealized example, let consider an ecosystem intrinsically poor in a key nutrient, such as nitrogen, as in many grasslands (e.g., ref. [START_REF] Abbadie | Lamto, Structure, Functioning and Dynamics of a Savanna Ecosystem[END_REF]). The conceptual model assumes part of the nitrogen (say, u) trapped in plant biomass and part (say, v) mobile in the soil. Plants capture N from v and release it through death and biomass decay. Local N concentrations of the u and v forms may vary in space, but the overall quantity remains constant over the system/landscape (an assumption that we will not ecologically discuss here).

From this framework of conservative interplay between mobile vs. fixed forms, the authors defined the "densityaggregation" model as belonging to a specific category (namely, "conservative reaction-diffusion") of the well-known reaction-diffusion models. u and v being the fixed and mobile fractions, the density-aggregation model stands as (after appropriate scaling):

[ 1]

As first specificity, both equations have the same reaction term, f(u,v) with converse signs, in accord with the conservation of u + v as a constant amount (τ). Therefore, f(u,v) includes both fluxes from v to u (demobilization, e.g., plant uptake) and fluxes from u to v (mobilization, e.g., via biomass decay). δ is the ratio of the u to v diffusion rates (i.e., diffusion ratio). In the field of reaction-diffusion systems, the authors placed the conservative system of type Eq. 1, in the class of "activator-depleted substrates" (u being called "activator" and v "substrate"). This can be illustrated in our idealized example of nitrogen-grass tussocks, simply saying that living plant biomass is autocatalytic ("activator") at least at low local density, depletes the N resource ("substrate") as it grows, and keeps it trapped until the plant tissues die and decay. The emergence of patterns requires as additional feedbacks that the activator (here plant biomass) reduces substrate release (e.g., through decreased tissue mortality) and/or amplifies depletion (e.g., through increased nutrient uptake).

Using qualitative analyses and complementary simulations of their general model, the authors established conditions for patterning from diffusive instability and drew characteristic properties of the patchy state, notably with respect to δ and τ as main driving parameters. Patches are here local maxima of the focal substance under its two forms. Without much loss of genericity, they specified (in Table 2) the form of a general model addressing examples of plant-nutrient and plant-sand dynamics (the "picky gazing" modeling is a little bit specific). In these two cases, the plant density is supposed proportional to u, which induces proportionality between plant and patch covers (e.g., in figures 4 and 5 in Siteur et al.). Moreover, patch-size distribution is expected to be log-normal while there is a slow, overall patch coarsening through time. This is due to the intrinsic instability of the smallest patches that progressively release the focal substance and thereby allow larger patches catching it. Interestingly, coarsening leaves unchanged the shape of the patch-size distribution (once rescaled by mean patch radius). Any local disturbance is speeding up coarsening since smallest patches never recover (figure 4). For this reason, the authors highlight the higher vulnerability of patchy systems originating in the density-dependent aggregation (DDA) process compared with outcomes of SDF. This lack of resilience at local scale is not expected for the entire system in which the focal substance is assumed conserved in coarser patches. The concerns of disturbance mitigation or ecosystem restoration are thus shifted toward the broadest scales. In the real world, this may raise questions about boundary conditions and minimum extents for system viability. [START_REF] Pringle | Spatial self-organization of ecosystems: Integrating multiple mechanisms of regular-pattern formation[END_REF], several patterning processes may simultaneously be acting. In the famous patterns of banded vegetation (tiger bush), some observers underlined the role of a "sandy veil" 21, which conservative dynamics may interact with the dissipative vegetation-water dynamics.

u t = f (u, v) + 2 u x 2 + 2 u y 2 v t = -f (u, v) + 2 v x 2 + 2 v y 2 A b) B
In real-world systems that are intrinsically heterogeneous and experiencing disturbances, tracking signatures of specific processes through indicators is naturally tricky, yet necessary as to nurture the fascinating dialogue between models and nature. A snapshot of a certain pattern (as in Fig. 1A here), whatever its impressiveness is rarely unequivocal about the underlying mechanisms. Field knowledge remains essential along with replicated observations that are fortunately increasingly available from aerial or space remote sensing. Moreover, observing large systems is often useful to access different stages of pattern development (e.g., Fig. 1B). The evaluation of well-identified model classes in light of rich multitemporal data sets is essential to allow ecological pattern science to progress to maturity and to provide insights into the resilience of ecosystems in the face of global change.

"The contribution of Siteur et al. is another example of the fruitfulness of introducing physical methods in ecology."

Fig. 1 .

 1 Fig. 1. Is this system dissipative with SDF or conservative with density-dependent aggregation? A-Top, Tussock grass pattern in the Andean Altiplano (Sajama National Park, Bolivia, see ref. 9 for context, photos: courtesy F. Anthelme). In this very case and apart from patch-size distribution, observations from a specific area (b-bottom) showing spiral patterns rather plead for localized structures and curvature instability (as in ref. 10) from SDF.

  The contribution of Siteur et al. (1), is another example of the fruitfulness of introducing physical methods in ecology. It installs DDA as a distinct patterning process, which may be of potential relevance to interpreting a vast array of ecological systems. Theoretical findings on patch-size distribution can provide clues that DDA is at play instead of other mechanisms, such as SDF or sole positive feedbacks. Coarsening may be more equivocal since it has already been encountered in some SDF models (see the "global competition model" in ref. 5 p.141, or ref. 8 relating to the Swift-Hohenberg equation). Moreover, in many real ecological systems (e.g., ref.

Downloaded from https://www.pnas.org by CIRAD on February 2, 2023 from IP address 194.199.236.40.


Marco Donnini

Jean-Luc Probst

Anne Probst

Francesco Frondini

Ivan Marchesini

Fausto Guzzetti

The contribution of weathering of the main Alpine rivers on the global carbon cycle

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The contribution of weathering of the main Alpine rivers on the global carbon cycle Marco Donnini (1,4), Jean-Luc Probst (2,3), Anne Probst (2,3), Francesco Frondini (1), Ivan Marchesini (4), and Fausto Guzzetti (4) The contribution of weathering of the main Alpine rivers on the global carbon cycle Abstract On geological time-scales the carbon fluxes from the solid Earth to the atmosphere mainly result from volcanism and metamorphic-decarbonation processes, whereas the carbon fluxes from atmosphere to solid Earth mainly depend on weathering of silicates and carbonates, biogenic precipitation and removal of CaCO 3 in the oceans and volcanic gases -seawater interactions. Quantifying each contribution is critical. In this work, we estimate the atmospheric CO 2 uptake by weathering in the Alps, using results of the study of the dissolved loads transported by 33 main Alpine rivers. The chemical composition of river water in unpolluted areas is a good indicator of surface weathering processes [START_REF] Garrels | Evolution of Sedimentary Rocks[END_REF][START_REF] Drever | The Geochemistry of Natural Waters[END_REF]Meybeck, 1984;[START_REF] Tardy | Cycle De l'Eau. Climats[END_REF][START_REF] Berner | The Global Water Cycle. Geochemistry and Environment[END_REF][START_REF] Probst | Carbon river fluxes and global weathering CO 2 consumption in the Congo and Amazon river basins[END_REF]. The dissolved load of streams originates from atmospheric input, pollution, evaporite dissolution, and weathering of carbonate and silicate rocks, and the application of mass balance calculations allows quantification of the different contributions. In this work, we applied the MEGA (Major Element Geochemical Approach) geochemical code (Amiotte Suchet, 1995; Amiotte Suchet and [START_REF] Amiotte-Suchet | Origins of dissolved inorganic carbon in the Garonne river waters: seasonal and interannual variations[END_REF] to the chemical compositions of the selected rivers in order to quantify the atmospheric CO 2 consumed by weathering in Alpine region. The drainage basins of the main Alpine rivers were sampled near the basin outlets during dry and flood seasons. The application of the MEGA geochemical consisted in several steps. First, we subtracted the rain contribution in river waters knowing the X/Cl (X = Na, K, Mg, Ca) ratios of the rain. Next, we considered that all (Na+K) came from silicate weathering. The average molar ratio Rsil = (Na+K)/(Ca+Mg) for rivers draining silicate terrains was estimated from unpolluted French stream waters draining small monolithological basins [START_REF] Meybeck | Composition Chimique Naturelle Des Ruisseaux Non Pollués En France[END_REF][START_REF] Meybeck | Global Chemical Weathering of Surficial Rocks Estimated from River Dissolved Load[END_REF]. For the purpose, we prepared a simplified geo-lithological map of Alps according to the lithological classification of [START_REF] Meybeck | Composition Chimique Naturelle Des Ruisseaux Non Pollués En France[END_REF][START_REF] Meybeck | Global Chemical Weathering of Surficial Rocks Estimated from River Dissolved Load[END_REF]. Then for each basin we computed Rsil weighted average considering the surface and the mean precipitation for the surface area of each lithology. Lastly, we estimated the (Ca+Mg) originating from carbonate weathering as the remaining cations after silicate correction. Depending on time-scales of the phenomena (shorter than about 1 million year i.e. correlated to the short term carbon cycle, or longer than about 1 million years i.e. correlated to the long-term carbon cycle), we considered different equations for the quantification of the atmospheric CO 2 consumed by weathering [START_REF] Huh | Estimation of Atmospheric CO 2 Uptake by Silicate Weathering in the Himalayas and the Tibetan Plateau: a Review of Existing Fluvial Geochemical Data[END_REF]. The results show the net predominance of carbonate weathering on fixing atmospheric CO 2 and that, considering the long-term carbon cycle, the amount of atmospheric CO 2 uptake by weathering is about one order of magnitude lower than considering the short-term carbon cycle. Moreover, considering the short-term carbon cycle, the mean CO 2 consumed by Alpine basins is of the same order of magnitude of the mean CO 2 consumed by weathering by the 60 largest rivers of the world estimated by [START_REF] Gaillardet | Global Silicate Weathering and CO 2 Consumption Rates Deduced from the Chemistry of Large Rivers[END_REF].

( 1 )

 1 Università degli Studi di Perugia. Dipartimento di Scienze della Terra. Perugia, Italy, (4) Consiglio Nazionale delle Ricerche, Istituto di Ricerca per la Protezione Idrogeologica, Perugia, Italy, (2) University of Toulouse ; INPT, UPS ; Laboratorie Ecologie Fonctionnelle et Environment (EcoLab), ENSAT, Castanet Tolosan, France, (3) Centre National de la Recherche Scientifique (CNRS), EcoLab, ENSAT Castanet Tolosan, France


MD Joanna Kedra

MD Raphaele Seror
email: raphaele.seror@aphp.fr

Philippe Dieudé

Arnaud Constantin
email: arnaud.constantin@univ-tlse3.fr

MD Eric Toussirot
email: etoussirot@chu-besancon.fr

MD Elias Kfoury
email: elias.kfoury@ch-brive.fr

MD Charles Masson
email: chmasson@chu-angers.fr

Divi Cornec
email: divi.cornec@chu-brest.fr

MD Jean-Jacques Dubost
email: jjdubost@chu-clermontferrand.fr

MD Laurent Marguerie
email: lmarguerie@hopale.com

Sebastien Ottaviani
email: sebastien.ottaviani@aphp.fr

MD Franck Grados
email: grados.franck@chu-amiens.fr

Rakiba Belkhir
email: rakiba.belkhir@aphp.fr

Olivier Fain

Bruno Fautrel

Peggy Philippe

Muriel Piperno

Bernard Combe

Olivier Lambotte

Christophe Richez

Jérémie Sellam

Thomas Sené

Guillaume Denis

Thierry Lequerre

Thierry Lazure

Xavier Mariette

Gaetane Nocturne

MD Philippe Dieude
email: philippe.dieude@aphp.fr

PhD Arnaud Md

MD Sébastien Ottaviani

Lack of association between the TNFAIP3 rs2230926 variant and rheumatoid arthritisassociated lymphoma

Abstract (215/250 words)

Objective: Given the association of TNFAIP3 rs2230926 variant with rheumatoid arthritis (RA) on the one hand, and Sjögren syndrome-related lymphoma on the other hand, the aim of this study was to investigate its contribution to RA-related lymphoma susceptibility.

Methods: In this multicentre case-control study, cases were patients with RA fulfilling ACR-EULAR 2010 criteria, who developed a B-cell non-Hodgkin lymphoma (NHL) or Hodgkin's lymphoma after the diagnosis of RA. Controls were patients with RA and without lymphoma matched on age. Blood or saliva samples were collected for genotyping.

Results: Among the 162 patients included in this study, DNA was available for 145 (89.5%) patients (38 cases and 107 controls), having no difference in demographic and disease characteristics compared to the whole study population. Overall, 100 patients (69.0%) were female, with a mean age of 51.9 years (SD=10.9) at RA diagnosis and a mean disease duration of 10.6 years (SD=5.0). Lymphomas were mostly diffuse large B-cell lymphoma (N=19, 50.0%), follicular lymphoma (N=6, 15.8%) and marginal zone lymphoma (N=4, 10.5%). The genotyping study revealed no significant association between the TNFAIP3 rs2230926 variant and the occurrence of lymphoma (p=0.30)

Conclusion:

We did not detect any association between the TNFAIP3 rs2230926 variant and RA-related lymphoma. This highlights the need for further studies to determine the specific mechanisms driving RA lymphomagenesis.

Key words: rheumatoid arthritis, lymphoma, TNFAIP3, A20, single nucleotide polymorphism

There is an increased risk of lymphoma in patients with auto-immune diseases (AID) especially rheumatoid arthritis (RA) and primary Sjögren Syndrome (pSS) with a respective relative risk of 2 and 15 [START_REF] Zintzaras | The risk of lymphoma development in autoimmune diseases: a meta-analysis[END_REF]. The mechanisms promoting lymphomagenesis in AID might be diverse [START_REF] Kedra | Inflammation-targeted therapies and cancer[END_REF]. Disease activity and chronic antigenic stimulation of autoimmune B cells play a key role in lymphomagenesis associated with pSS [START_REF] Nocturne | Lymphomas complicating primary Sjögren's syndrome: from autoimmunity to lymphoma[END_REF]. In this context, a complete control of the NF-kB pathway activation is required to avoid B-cell clonal proliferation. TNFAIP3 encodes for the protein A20, a key gatekeeper of NF-kB activation and the TNFAIP3 rs2230926 variant was previously found to be associated with lymphoma in pSS [START_REF] Nocturne | Germline and somatic genetic variations of TNFAIP3 in lymphoma complicating primary Sjogren's syndrome[END_REF][START_REF] Nocturne | Germline variation of TNFAIP3 in primary Sjögren's syndrome-associated lymphoma[END_REF]. This variant is responsible for a slight decreased control of NF-kB activation. In patients with RA, mechanisms leading to lymphoma occurrence, mainly diffuse large B-cell lymphoma (DLBCL), is less understood. From an epidemiologic point of view, it has been demonstrated that high cumulative RA disease activity increase odds for lymphoma [START_REF] Baecklund | Association of chronic inflammation, not its treatment, with increased lymphoma risk in rheumatoid arthritis[END_REF]. Given the association of TNFAIP3 rs2230926 with the overall RA [START_REF] Musone | Sequencing of TNFAIP3 and association of variants with multiple autoimmune diseases[END_REF], and pSS-related lymphoma [START_REF] Nocturne | Germline and somatic genetic variations of TNFAIP3 in lymphoma complicating primary Sjogren's syndrome[END_REF][START_REF] Nocturne | Germline variation of TNFAIP3 in primary Sjögren's syndrome-associated lymphoma[END_REF], we decided to investigate its contribution to RA-related lymphoma susceptibility.

To this end, a multicentre case-control study was performed. Cases were patients with RA fulfilling ACR-EULAR 2010 criteria, who developed a B-cell NHL or Hodgkin's lymphoma after the diagnosis of RA. Controls were patients with RA and without lymphoma matched on age. Fifty-four cases of RA-related lymphoma were included and matched to 108 controls (1:2) [START_REF] Kedra | Lymphoma complicating rheumatoid arthritis: results from a French case-control study[END_REF]. DNA was available for 145/162 (89.5%) patients (38 cases and 107 controls), having no difference in demographic and disease characteristics compared to the whole study population (table 1). Overall, 101 patients (69.2%) were female, with a mean age of 51.9 years (SD=10.9) at RA diagnosis and a mean disease duration of 10.5 years (SD=4.9). Lymphomas were mostly DLBCL (N=19, 50.0%), follicular lymphoma (N=6, 15.8%) and marginal zone lymphoma (N=4, 10.5%). The genotyping study revealed no significant association between the TNFAIP3 rs2230926 variant and the occurrence of lymphoma (p=0.30) (table 2).

Even if this study could suffer from a relatively low power of detection (estimated power of 80%), this lack of association suggests distinct mechanisms of lymphomagenesis in RA and pSS. Lymphoma histology and localization differ in pSS and RA associated lymphoma: MALT lymphoma of salivary glands in pSS versus nodal DLBCL in RA. RF positivity and/or cryoglobulinemia are associated with lymphoma occurrence in pSS but not in RA [START_REF] Nocturne | Lymphomas complicating primary Sjögren's syndrome: from autoimmunity to lymphoma[END_REF]. In RA, instead of a direct transformation of autoimmune B cells, lymphoma might derive from associated subsets of B cells such as memory B cells. Interestingly, memory B cells involved in RA are increased in longstanding disease [START_REF] Fedele | Memory B cell subsets and plasmablasts are lower in early than in long-standing rheumatoid arthritis[END_REF]. Recently, TBL1XR1 rare variants were found to induce ABC DLBCL by promoting aberrant memory B cells, which are prone to avoid plasmocytic differentiation on recall and to undergo systemic dissemination [START_REF] Venturutti | TBL1XR1 Mutations Drive Extranodal Lymphoma by Inducing a Pro-tumorigenic Memory Fate[END_REF].

In aggregates, we did not detect any association between the TNFAIP3 rs2230926 variant and RA-related lymphoma. This highlights the need for further studies to determine the specific mechanisms driving RA lymphomagenesis. Footnote: IQR: interquartile range; ACPA: anti-citrullinated peptide antibodies; DAS28: disease activity score in 28 joints; RA: rheumatoid arthritis; RF: rheumatoid factor; the time of matching corresponded to the diagnosis of lymphoma for cases and the 10-year ESPOIR visit for controls

Table 1 :

 1 characteristics of the genotyping study population

		Whole genotyping
		study population	Cases (N=38)	Controls (N=107)
		(N=145)
	Male gender, N (%)	45 (31.0)	20 (52.6)	25 (23.4)
	Age at RA diagnosis, mean (SD)	51.9 (10.9)	49.8 (12.6)	52.7 (10. 2)
	Age at the time of matching, mean (SD)	62.0 (10.1)	61.2 (10.2)	62.2 (10.1)
	ACPA positivity, N (%)	106 (73.1)	35 (92.1)	71 (66.4)
	RF positivity, N (%)	111 (76.6)	35 (92.1)	76 (71.0)
	Erosions on X-rays at
	the time of matching, N	55 (37.9)	30 (78.9)	25 (23.4)
	(%)
	DAS28 at the time of matching, mean (SD)	2.9 (1.5)	3.9 (1.7)	2.6 (1.4)
	RA duration at the time of matching, mean (SD)	10.6 (5.0)	11.9 (9.7)	10.0 (0.3)
	Year of RA diagnosis, median (IQR)	2002 (1978-2016)	2000 (1978-2016)	2004 (2002-2005)
	Year of lymphoma
	diagnosis/10-year ESPOIR visit, median	2013 (1996-2018)	2012 (1996-2018)	2014 (2012-2015)
	(IQR)

Table 2 :

 2 Testing for association of TNFAIP3 rs2230926 with RA-related lymphoma

	TNFAIP3	Cases (N=38)	Controls (N=107)	OR (95%	p-value*
	rs2230926			CI)
	MAF	2.6%	5.6%	0.45 (0.10-	0.30
				2.08)
	TT	36 (94.7%)	96 (89.7%)	0.49 (0.10-	0.35
				2.30)
	TG	2 (5.3%)	10 (9.3%)	0.53 (0.11-	0.42
				2.55)
	GG	0 (0.0%)	1 (0.9%)	-	0.54
	MAF: minor allele frequency
	*chi-square test

Acknowledgements:

The authors wish to acknowledge the French Society of Rheumatology (Société Française de Rhumatologie -SFR) for the grant provided to support the genetic analysis.

Statements and declarations:

ESPOIR cohort: An unrestricted grant from Merck Sharp and Dohme (MSD) was allocated for the first 5 years. Two additional grants from INSERM were obtained to support part of the biological database. The French Society of Rheumatology, Abbvie, Pfizer, Lilly, and more recently Fresenius and Galapagos also supported the ESPOIR cohort study. We also thank Nathalie Rincheval for expert monitoring and data management and all the investigators who recruited and followed the patients (F. BERENBAUM, Paris-Saint Antoine, MC. BOISSIER, Paris-Bobigny, A .CANTAGREL, A CONSTANTIN, Toulouse, B. COMBE , Montpellier, M.DOUGADOS, Paris-Cochin, P BOUMIER, Amiens, B. FAUTREL, Paris-La Pitié, RM. FLIPO, Lille , Ph. GOUPILLE, Tours, F. LIOTE, Paris-Lariboisière, O VITTECOQ, Rouen, X MARIETTE, Paris Bicetre, Ph DIEUDE, Paris Bichat, A.SARAUX, Brest, Th SCHAEVERBEKE, Bordeaux, J. SIBILIA, Strasbourg). One biological resources centre (Paris-Bichat, S Tubiana) was in charge of centralising and managing biological data collection.

Funding details:

The French Society of Rheumatology (Société Française de Rhumatologie -SFR) provided a grant to support the genetic analysis.

Competing interests:

The authors have no disclosure to declare in relationship to the present study.


Abdelaati Daouia

Gilles Stupfler

Antoine Usseglio-Carleve

Extreme value modelling of SARS-CoV-2 community transmission using discrete Generalised Pareto distributions

Keywords: COVID-19, Superspreading, Cluster size, Secondary cases, Extreme value theory, Discrete extremes

Extreme value modelling of SARS-CoV-2 community transmission using discrete Generalised Pareto distributions

Abdelaati Daouia, Gilles Stupfler, Antoine Usseglio-Carleve

Introduction

Superspreading events (SSEs) have been recognised as a significant source of disease transmission for respiratory coronaviruses such as SARS-CoV and SARS-CoV-2 [START_REF] Lloyd-Smith | Superspreading and the effect of individual variation on disease emergence[END_REF][START_REF] Adam | Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong[END_REF]. SSEs may be defined as outbreaks in which a given individual (the index case) infects a number of people (secondary cases) well above a certain measure, such as the average or median number of infections. The number of secondary cases resulting directly from an index case can be viewed as a random variable, say Z, defining the so-called offspring distribution. For both coronaviruses, events having triggered more than 6 secondary cases have been suggested to constitute SSEs [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF].

Data on such SSEs that was curated and reported in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] in the early stages of the COVID-19 pandemic is necessarily scarce: it consists mainly of 15 SSEs associated with SARS-CoV and 45 SSEs associated with SARS-CoV-2, each represented by a number of secondary cases Z i resulting from a single given index case in Europe, Asia or North America. The natural framework for the analysis of SSEs, and more generally of atypical observations far away from the mean, is extreme value theory.

Following this framework, it was argued in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] that SSEs are fat-tailed, although this was done by pooling the 60 available SSEs from SARS-CoV and SARS-CoV-2. A careful investigation of these SARS-CoV and SARS-CoV-2 datasets reveals that the two largest observations in the pooled data are SARS-CoV SSEs; given the small sample size, one may wonder whether the reported estimate of tail heaviness is representative of the tail behaviour of SARS-CoV-2 SSEs. This constitutes the motivation for this work, whose overarching goals are to show how to conduct a statistically rigorous extreme value analysis of community transmission parameters, and to carry out such an analysis in the example of SARS-CoV-2. By focusing directly on the raw SARS-CoV-2 data considered in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF], we provide evidence of a lighter upper tail for SSEs with significantly less tail exposure than predicted in their study. We arrive at the same conclusion by making use of a more recent and much larger publicly available surveillance and contact-tracing database containing the number of secondary cases Z i for 88,527 index cases in the Indian states of Andhra Pradesh and Tamil Nadu [START_REF] Laxminarayan | Epidemiology and transmission dynamics of COVID-19 in two Indian states[END_REF]. We also analyse two other South Korean contact-tracing datasets, one collected in the first half of 2020 [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF], the other during the summer of 2021 when the Delta variant of SARS-CoV-2 was responsible for the majority of positive cases [START_REF] Ryu | Serial interval and transmission dynamics during SARS-CoV-2 Delta variant predominance, South Korea[END_REF]. The fat-tailedness of the secondary cases distribution is found to be rather clear in the 2021 sample of data, while the analysis of the 2020 data is less conclusive. In all these samples of data we find point estimates of the extreme value index suggesting that the secondary cases distribution has a finite third moment, which stands in contrast with the earlier finding of [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] of a distribution with an infinite variance.

In addition to that, we consider the broader problem of large community transmission, as it represents the other fundamental source of pandemic risk. Large infection clusters, along with SSEs, have been argued to play an important role in the transmission of SARS-CoV-2 [START_REF] Adam | Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong[END_REF]. In a similar spirit to [START_REF] Adam | Clustering and superspreading potential of SARS-CoV-2 infections in Hong Kong[END_REF], we define a cluster of SARS-CoV-2 cases in our analysis as a local outbreak involving a minimum of two cases, including confirmed close contacts with epidemiological linkage over a limited period of time. We consider two databases constructed from government reports [6, 7, 8, 9] and media sources [START_REF] Swinkels | SARS-CoV-2 Superspreading Events Database[END_REF], comprising 15 samples of SARS-CoV-2 cluster sizes recorded in 11 countries and 4 US states. Our results show that 13 of these 15 countries and states have fat-tailed cluster size distributions, thus facilitating the process of inferring their risk category in terms of large community transmission. This allows us to better understand the drivers of superspreading and cluster formation in the ongoing COVID-19 pandemic. The recent theory of discrete extremes [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF][START_REF] Prieto | Modelling road accident blackspots data with the discrete generalized Pareto distribution[END_REF][START_REF] Hitz | Discrete Extremes[END_REF][START_REF] Ranjbar | Modelling the extremes of seasonal viruses and hospital congestion: the example of flu in a Swiss hospital[END_REF] is our basic tool to address the highly discrete nature of SARS-CoV-2 secondary transmission data and cluster sizes. Its use constitutes our main statistical contribution to the study of the transmission of the SARS-CoV-2 virus. As we illustrate throughout the paper, estimating and inferring the extreme value index and extreme percentiles of the underlying discrete distributions with this methodology is much easier and accurate than with classical extreme value methods such as the Hill and Generalised Pareto maximum likelihood estimators, which heavily rely on the continuous data assumption.

The structure of the paper is as follows. We first describe the methods employed throughout our study, including the discrete Generalised Pareto Distribution fitted to exceedances over a high threshold by means of the maximum likelihood estimator.

We then analyse our datasets, first on SARS-CoV-2 secondary case numbers and then on cluster sizes, using these methods. The final section gathers and contrasts these findings and concludes with additional comments about the scope, limitations and robustness of our results, as well as ideas for further work.

Methods

We use several methods from extreme value theory, which constitutes the correct mathematical framework for the analysis of high observations from a random phenomenon [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]. We are particularly interested in methods that can describe so-called fat-tailed random variables, which infrequently but regularly generate very high values and therefore appear to be relevant in the analysis of SARS-CoV-2 transmission. A random variable X is fat-tailed if and only if its distribution function P(X ≤ x) can be, for large x, expressed as P(X ≤ x) = 1 -x -1/ξ (x), where satisfies (tx)/ (t) → 1 as t → ∞ for any positive real number x. Informally, the tail behaviour of X is controlled by the extreme value index ξ > 0, which must be estimated to get a precise understanding of tail heaviness. A standard estimator in this context is the Hill estimator [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF]. For a dataset Z 1 , . . . , Z n , the Hill estimator at threshold u is defined as

ξH u = 1 n i=1 1 {Z i >u} n i=1 log Z i u 1 {Z i >u} .
It is of course crucial, before using the Hill estimator, to ascertain whether the distribution of the data points indeed has a heavy tail. A common diagnostic method is the mean excess plot, which estimates the values of the mean excess function

E(u) = E [Z -u|Z > u] as function of u. A natural estimate of E(u) is given, for each threshold u, by its empirical counterpart Ê(u) = n i=1 Z i 1 {Z i >u} n i=1 1 {Z i >u} -u.
A fat-tailed distribution will typically have mean excess plots exhibiting a linear upward drift for large values of u, whose slope is a consistent estimate of ξ/(1 -ξ) when ξ < 1, see for example Section 1.2.2 pp.14-19 and p.152 in [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF]. In the case ξ ≥ 1, Theorem 3.4 and Remark 3.5 in [START_REF] Ghosh | A discussion on mean excess plots[END_REF] show that the mean excess plot converges in a suitable sense to a random curve, which in the log-log scale is a straight line with slope 1/ξ and random intercept term constructed upon a stable random variable with index 1/ξ.

It has, however, been observed in the extreme value literature [START_REF] Ghosh | A discussion on mean excess plots[END_REF] that the mean excess function very often exhibits a non-linear behaviour at the right end of the mean excess plot, due to very high variability of the estimate of E(u) when u is close to the highest Z i . As a consequence, good statistical practice recommends to confirm a diagnostic of a heavy tail using other extreme value tools. One such general approach, which does not presuppose that the data is fat-tailed, consists in using the Generalised Pareto maximum likelihood estimator applied to the excesses

Z i -u.
Recall that the Generalised Pareto distribution, with shape parameter ξ and scale parameter σ, has probability density function

h ξ,σ (x) = 1 σ 1 + ξ x σ -1/ξ-1 1 {x>0, 1+ξx/σ>0} .
The Generalised Pareto maximum likelihood estimator is then defined as, according to Section 5.3.2 in [START_REF] Beirlant | Statistics of Extremes: Theory and Applications[END_REF] and [START_REF] Zhou | The extent of the maximum likelihood estimator for the extreme value index[END_REF]:

ξGP u , σGP u = arg min ξ>-1, σ>0 n i=1 log h ξ,σ (Z i -u) = arg min ξ>-1, σ>0 n i=1 -log σ - 1 ξ + 1 log 1 + ξ Z i -u σ 1 {Z i >u, 1+ξ(Z i -u)/σ>0} .
The Generalised Pareto maximum likelihood estimators are valid even when the underlying distribution is not fat-tailed, which has made them very popular in the natural sciences [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF].

However, both the Hill and Generalised Pareto estimators of ξ suffer from jagged sample paths when the data points Z i feature a substantial number of ties, that is, they come from a distribution with a high degree of discreteness. This behaviour makes it extremely difficult to choose an accurate estimate of ξ, which renders the two methods highly unsatisfactory. The essential reason behind this phenomenon is that both estimators are built under the -generally incorrect -assumption that the data points come from a pure (Generalised) Pareto distribution, which is continuous, and as such, they cannot be expected to handle a substantial degree of discreteness. We exemplify this phenomenon in Fig. 1: notice, in the top panels, the stark difference in stability and smoothness of sample paths between a plot of the Hill estimator as a function of the threshold value (henceforth referred to as a Hill plot) for continuous data Z i and its counterpart for data rounded to the nearest integer up. Crucially in applied setups, the asymptotic Gaussian confidence intervals constructed by approximating the distribution of nP(X > u n )( ξH un -ξ) by a Gaussian distribution with expectation 0 and variance ξ 2 , which is valid when u n → ∞ satisfies reasonable conditions [START_REF] Goldie | Slow variation with remainder: Theory and applications[END_REF], are highly unstable when the data features a large number of ties, thus making inference using the Hill estimator unadvisable.

The bottom panels further show the impact of these data ties: the Hill estimator for Hill plots (solid lines) and corresponding 90% Gaussian asymptotic confidence intervals (dashed lines) as functions of the threshold value u, for n = 10,000 simulated data points Z i from the Burr distribution with probability density function f (x) = ξ -1 x -ρ/ξ-1 (1 + x -ρ/ξ ) 1/ρ-1 (for x > 0) with ξ = 1/2 and ρ = -1 in the left panel, and for the data Z i (i.e. the smallest integer larger than or equal to Z i ) in the right panel. Bottom panels: Averaged Hill plots when this experiment is repeated N = 1,000 times.

discrete data tends to be strongly biased and much more so than the Hill estimator for continuous data.

An alternative option properly taking the discreteness of the data into account is to employ discrete models to construct an estimator of the extreme value index.

This was pursued by [START_REF] Shimura | Discretization of distributions in the maximum domain of attraction[END_REF][START_REF] Hitz | Discrete Extremes[END_REF], which used so-called D-GPD (for Discrete-Generalised Pareto Distribution) models, first employed by [START_REF] Prieto | Modelling road accident blackspots data with the discrete generalized Pareto distribution[END_REF] to model road accidents and more recently in [START_REF] Ranjbar | Modelling the extremes of seasonal viruses and hospital congestion: the example of flu in a Swiss hospital[END_REF] to model hospital congestion. The D-GPD, whose probability mass function is When ξ = 0, the convention we adopt is that (1 + ξz) -1/ξ = exp(-z), for any z ∈ R. These maximum likelihood estimators of the extreme value index ξ and scale parameter σ of the D-GPD model are readily obtained through the R maximisation routine optim. Using the classical theory of maximum likelihood estimators, confidence intervals for ξ may be derived from ξu by estimating the total Fisher information matrix I(ξ, σ) using a finite difference method and then deducing the following 100α%-confidence interval for ξ:

p ξ,σ (x) = 1 + ξ x σ -1/ξ -1 + ξ x + 1 σ -1/ξ for x = 0,
ξu + Î(ξ, σ) -1 1,1 Φ -1 1 -α 2 , ξu + Î(ξ, σ) -1 1,1 Φ -1 1 + α 2 ,
where Φ denotes the standard normal distribution function and Φ -1 its inverse (quantile function). Modelling Z -u conditional on Z ≥ u by a D-GPD distribution with parameter estimates ( ξu , σu ) suggests the following estimate of the 100αth percentile of Z adapted from [12, Formula (5) p.41]:

qα = σu ξu n(1 -α) n i=1 1 {Z i ≥u} -ξu -1 + u -1 ,
for α ∈ (0, 1) large enough. Here, • denotes the ceiling function, that is, x denotes the smallest integer larger than or equal to x. Estimating this quantile by plugging in the aforementioned estimates of ξ and σ makes it possible to infer extreme quantile levels and therefore get precise information on the tail behaviour of a distribution with a large degree of discreteness. For each of the extreme value estimators we have introduced (Hill estimator, GPD and D-GPD maximum likelihood estimators), a common practice for selecting a suitable pointwise estimate of ξ, colloquially referred to as "eyeballing", is to pick out a sufficiently high threshold u corresponding to a stable region of the plot [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF]. We shall indeed also adopt this practice and will clearly indicate selected thresholds or threshold regions in our analyses.

For comparison purposes, we will contrast the resulting extreme quantile estimates with those provided by the (conditioned) negative binomial distribution. Recall that the probability mass function of the negative binomial distribution (with parameters r > 0 and p ∈ (0, 1)) conditional on Z > u, is given by

P p,r,u (Z = k) = Γ(k+r) k! Γ(r) p r (1 -p) k 1 -u i=0 (Γ(i + r)/(i! Γ(r)))p r (1 -p) i , for all k > u.
Here Γ denotes Euler's Gamma function. With a dataset z 1 , . . . , z n , the parameter estimators are therefore obtained as the maximum log-likelihood solution arg max

(p,r)∈(0,1)×(0,∞) n i=1 log P p,r,u (Z = z i ).
Ever since the seminal work of [START_REF] Lloyd-Smith | Superspreading and the effect of individual variation on disease emergence[END_REF], the negative binomial distribution has been widely used to describe the number of secondary cases resulting from an index case of SARS-CoV. As suggested in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF][START_REF] Kremer | Quantifying superspreading for COVID-19 using Poisson mixture distributions[END_REF], this model has exponentially decreasing probability mass functions and thus cannot be expected to accurately represent tail heaviness in SARS-CoV-2 transmission data. We provide below further evidence for this claim, and for the suitability of D-GPD maximum likelihood estimates in the context of discrete data, through several datasets gathering numbers of SARS-CoV-2 secondary cases and cluster sizes in different settings.

Data and results

Analysis of secondary case data. Our first two datasets were reported in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF].

They consist of 15 SSEs associated with SARS-CoV (Dataset S1) and 45 SSEs associated with SARS-CoV-2 (Dataset S2), each resulting in more than 6 secondary cases, along with month of occurrence and location of the superspreading event, and its setting. We refer to [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] for further details about the construction of these datasets. Pooling the 15 SSEs associated with SARS-CoV and 45 SSEs associated with SARS-CoV-2 into a single sample and making use of a Generalised Pareto approximation, [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] has suggested that the distribution of the number of secondary cases Z belongs to the Fréchet maximum domain of attraction [START_REF] Cirillo | Tail risk of contagious diseases[END_REF], that is, the set of Pareto-type distributions, with extreme value index ξ between 0.5 and 1 (the estimate provided in [3, Fig. 1 E] is ξ ≈ 0.6). The index ξ tunes the tail heaviness of the distribution, with higher positive values indicating a heavier upper tail: moments of order higher than or equal to 1/ξ do not exist. An estimate of ξ around 0.6 means that the second moment of Z does not exist, reflecting the outsized contribution of SSEs to overall transmission. Most importantly perhaps, these findings on the tail heaviness of Z invalidate the conventional assumption that Z follows a negative binomial distribution for either coronavirus, whereas this assumption was widely adopted in the literature on disease transmission ever since the influential work [START_REF] Lloyd-Smith | Superspreading and the effect of individual variation on disease emergence[END_REF] on SARS-CoV, and it is still widely employed for SARS-CoV-2, see e.g. [START_REF] Ryu | Serial interval and transmission dynamics during SARS-CoV-2 Delta variant predominance, South Korea[END_REF][START_REF] Islam | COVID-19 and climatic factors: A global analysis[END_REF][START_REF] Hwang | Transmission dynamics of the Delta variant of SARS-CoV-2 infections in South Korea[END_REF].

Based on our statistical analysis of these datasets, summarised in Fig. 2, one may however argue that the method of [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] is inappropriate for examining the tail behaviour of their particular 60 SSEs. The sparsity of data on SSEs is addressed by combining the 15 and 45 observations associated with SARS-CoV and SARS-CoV-2 into a single sample, whereas the two datasets correspond to completely different distributions (Fig. 2 (a)) and should not be pooled accordingly. This is apparent from either a Kolmogorov-Smirnov test, with p-value 0.015, or the more common approach making the questionable assumption that Z follows a negative binomial distribution. The conditional (given Z > 6) negative binomial fit of the probability mass function to the Z i (by construction larger than 6), calculated as described in the last paragraph of the Methods section (Fig. 2 (b)), already suggests that the upper tail of Z for SARS-CoV appreciably dominates that for SARS-CoV-2.

In other words, even a naive analysis of the SSE distributions, using the classical negative binomial distribution and not accounting for the heavy tail in the data, indicates that the SSEs for SARS-CoV and those for SARS-CoV-2 exhibit different statistical behaviour. This is confirmed by an analysis of the data properly taking into account its extremes (Fig. 2 (c)): the ξ estimates obtained from the Hill estimator in the special case of SARS-CoV-2 vary between 0.35 and 0.45, and as such differ substantially from the various competing estimates found to vary between 0.5 and 1 in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF]. Even the 90% confidence intervals of ξ for SARS-CoV-2 (dashed red lines in Fig. 2 (c Yet, the low sample size of this SSE dataset puts a question mark over the quality of the statistical analysis. Trustworthy extreme value inference may require a larger sample size, of the order of at least several thousands. This is why we also analysed a much larger Indian secondary case dataset of size n = 88,527 (Database S3). This comprehensive surveillance and contact-tracing database was collected in 2020 by the public health authorities of the two Indian states of Andhra Pradesh and Tamil Nadu, whose residents total about 10% of India's population. It was studied for instance in [START_REF] Laxminarayan | Epidemiology and transmission dynamics of COVID-19 in two Indian states[END_REF] and [START_REF] Kremer | Quantifying superspreading for COVID-19 using Poisson mixture distributions[END_REF], and we refer to the latter for more information about the database's construction and contents. Results are reported in Fig. 3. Although the barplot of this data (Fig. 3 (a)) gives evidence of a considerable right skewness and its summary extreme value analysis (Fig. 3 (b)) suggests a heavy right tail, it should be noted that since the Z i range from 0 to 39 with a sample size of 88,527, the data is necessarily highly discrete with a large number of tied observations (see Table 1). Ignoring the discrete nature of the Z i by modelling their tail behaviour with the (Generalised) Pareto distribution is inappropriate as this typically results in unreliable extreme value index estimates and confidence intervals [START_REF] Hitz | Discrete Extremes[END_REF]. This becomes obvious here by superimposing both the classical Hill and continuous Generalised Pareto maximum likelihood estimators of the extreme value index, as functions of a varying threshold u in Fig. 3 (c). Clearly, both plots are so volatile and jagged that it is hard to identify any stable region and therefore a reasonable point estimate of ξ cannot easily be determined. Using the D-GPD distribution to fit exceedances Z i -u above the threshold u (rather than trying to fit the whole of the distribution, as [START_REF] Kremer | Quantifying superspreading for COVID-19 using Poisson mixture distributions[END_REF] did using a discrete Pareto distribution) results in a much smoother and stable fit (Fig. The first dataset, which consists of n = 5,165 numbers of SARS-CoV-2 secondary cases Z i , was analysed in [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF]. See Table 2.

We revisit the estimation of, and inference about, the underlying extreme value index by comparing the D-GPD estimates with the classical GPD and Hill estimates.

Results are displayed in Fig. 4. A least squares fit to the first part of the mean excess plot (Fig. 4 (b)) suggests a linearly increasing fit to the mean excess function with a slope of around 0.85, but this ignores the flat or even slightly linearly decreasing right-hand part of the data cloud. This throws the assumption that the offspring distribution is fat-tailed in doubt, although the barplot of the data (Fig. 4 (a)) would tentatively back the heavy tail assumption. The Hill estimator, which presupposes that the data is fat-tailed and graphed as a black line in Fig. 4 (c), does not exhibit any stable region which would allow to produce a reasonable point estimate.

In such scenarios, best practice in extreme value theory requires calculating alternative extreme value estimators whose consistency does not rest upon the heavy tail assumption (unlike the Hill estimator), such as the general GPD and D-GPD estimators. These are also represented in Fig. 4 (c). Clearly, the paths of these two estimates follow a similar trajectory which is very different from that of the Hill plot. They point towards substantially lower estimates of ξ, and even though the estimates are overall larger than 0, the validity of the heavy tail assumption ξ > 0

is not obvious for this dataset. Fig. 4 (d) further supports this observation: in the (somewhat) stable region around the threshold u = 10, the 90% confidence interval produced through maximum likelihood theory contains the value 0. Our conclusion from the analysis of this dataset is that the distribution of the number of secondary cases is either fat-tailed but with a low extreme value index, or perhaps even has an exponential-type tail. As a consequence, our finding is qualitatively different from that of [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF], since we do not obtain ξ estimates similar to those found by merging Datasets S1 and S2.

The second South Korean contact-tracing dataset comprises n = 33,903 SARS-CoV-2 numbers of secondary cases Z i (Database S5) detected between 25th July 2021 and 15th August 2021. It was initially explored in [START_REF] Ryu | Serial interval and transmission dynamics during SARS-CoV-2 Delta variant predominance, South Korea[END_REF], where it was highlighted that the Delta variant accounted for the majority of those cases. We therefore investigate this dataset to ascertain whether the tail behaviour of SSEs is substantially different for the Delta variant. The data is presented in Table 3 below. The results we obtain for this dataset are displayed in Fig. 5. The barplot of the data in (d) Discrete GPD maximum likelihood estimates of ξ (solid red) and their associated 90% confidence intervals (dashed red). In panels (c) and (d), the averaged discrete GPD estimate ξ = 0.209 over the stable region u ∈ [1, 10] is indicated with the horizontal red line.

Analysis of cluster size data. We broaden our analysis by examining whether SARS-CoV-2 cluster sizes are fat-tailed. We consider a database of 15 samples of cluster sizes recorded in 11 countries and 4 US states. We define a cluster as a local outbreak involving a minimum of two cases, including confirmed close contacts with epidemiological linkage observed up to extinction of the outbreak. This differs from the number of secondary cases linked to a single, given index case in an SSE, since the cluster size is now the total number of infected people over the duration of the outbreak. The number of reported clusters per country or state varies from 29 (France) to 4,769 (Colorado, USA). The database is constructed from government reports [6, 7, 8, 9] (Database S6) and media sources [START_REF] Swinkels | SARS-CoV-2 Superspreading Events Database[END_REF] (Database S7). The median cluster sizes were 5 (Database S6) and 33 (Database S7), and the largest clusters had sizes 1,761 (Database S6, in a Colorado prison) and 7,000 (Database S7, in an Italian football stadium). We denote by Y i the number of SARS-CoV-2 cases in cluster i.

The ξ estimates from each sample of cluster sizes allow to infer the risk category of the corresponding country/state in terms of local community transmission. 4, where 13 out of the 15 countries or states appear to have fat-tailed cluster size distributions (confirmed at the 90% confidence level except for China). We note that there is strong variation in point estimates of ξ across countries and states. The low sample sizes of the data available in each case (except for the two US states of Colorado and Oregon) certainly play an important role in that variation. Heterogeneity in population density and healthcare policies may also be substantial factors, although this would have to be cross-checked using complete demographic and public health data. The analysis for California and the United Kingdom was inconclusive. For the California dataset, this is possibly due to a strong degree of heterogeneity (see the histogram in the bottom left panel of Fig. 7). A stratified study of the Californian data might be more conclusive. For the UK dataset, the fact that the sample is so small (26 clusters) in a country with a highly developed healthcare and contact tracing system is suspicious and may suggest reporting issues.

Using the D-GPD model, one can gain further insight into large cluster sizes by providing extrapolated estimates of extreme percentiles q α potentially beyond the sample maximum, through the estimate qα described in the Methods section.

Estimated 95th and 99th percentiles are given in Table 4. One may also match 6, with the results obtained from the data whose sources were the media (Database S7). The top 9 plots refer to those countries and states for which the extreme value analysis was conclusive. The bottom 2 plots refer to those for which the extreme value analysis was inconclusive.

Discussion

In summary, we have investigated four datasets of secondary case numbers Z i for SARS-CoV-2 as a way to estimate and infer the extreme value index of the related underlying offspring distribution. Motivated by the highly discrete nature of such data, we used the Discrete GPD (D-GPD) maximum likelihood estimation method which produces smoother and more stable plots of the associated D-GPD estimator than the classical continuous GPD and Hill estimators. We first provided evidence that the small SSE dataset (Dataset S2) compiled by [START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF] We broaden our analysis by providing evidence that SARS-CoV-2 cluster sizes are typically fat-tailed, based on 15 samples from 11 countries and 4 US states.

We infer the risk exposure and risk category of each country and state by making use of D-GPD maximum likelihood estimates of both the extreme value index and extreme percentiles, along with their associated confidence intervals. For the sake of simplicity, we used a straightforward threshold selection rule, which is to spot a stability region in the estimates (as a function of the threshold value) and choose an estimate whose value is representative of those reached in this region. This practice, colloquially known as "eyeballing", is standard in applied extreme value analysis: A limitation of our study lies in the quality of the data, as it is not obvious whether all SSEs or clusters over a given time period were available, or whether cluster sizes were correctly recorded. To check robustness against missing data, we have reproduced part of our analysis of cluster data by removing 10% of observations at random in each sample containing at least 100 data points, and replicating this experiment 10,000 times. Robustness against poor recording was checked by multiplying each observation Y i by an independent normal variate W i having mean µ = 1 and standard deviation σ = 0.05, and then reproducing our analysis of cluster data on the Y i = W i Y i , this experiment being again replicated 10,000 times. There is indeed some variation in the resulting estimates of ξ (Figs. 8 and9), but this does not affect our conclusion on the fat-tailed behaviour of the data, except in rare situations when almost all the large values in the data go missing. This highlights the importance of accurate data reporting as a prerequisite to such analyses. A further limitation lies in the assumption of independent data that is implicitly made in order to derive confidence intervals for extreme value parameters, even though the data are implicitly time series. Handling serial dependence in the current setting of discrete epidemiological data is obviously an interesting but very difficult question,

involving the hitherto open problem of extreme value dependence in discrete time series, which deserves a study of its own.

It should be noted that, in classical epidemiological models, accurate estimation of the basic reproduction number R 0 is of crucial importance as it informs the extent of restrictions on social interactions and other control measures that should be imposed to terminate the spread of an epidemic. The range of R 0 for SARS-CoV-2 has been revised in [START_REF] Kochańczyk | Super-spreading events initiated the exponential growth phase of COVID-19 with R 0 higher than initially estimated[END_REF] to 4.7-11.4, which is considerably higher than most early estimates. This might explain why moderate restrictions that were implemented in some nations, e.g. France, Italy, Spain, the UK, Australia and New Zealand, turned out to be insufficient and replaced by nationwide or statewide lockdowns and/or border closures. It should be clear that our results are, by construction, robust to misspecified estimates of the expected number of secondary cases R 0 since they solely rely on extreme values of numbers of secondary cases.

Our approach can be viewed as a proof of concept that transmission data from a respiratory disease should not be pooled with data from a similar disease, since similar R 0 numbers or parameters of average transmission do not, in general, induce similar parameters of large community transmission. As such, preparing proactive control measures actually requires a fine assessment of how unequal the distributions of SSEs associated with different SARS-CoV-2 variants are. [START_REF] Liu | The reproductive number of the Delta variant of SARS-CoV-2 is far higher compared to the ancestral SARS-CoV-2 virus[END_REF] conclude that the reproductive number of the Delta variant is far higher than that of the historical SARS-CoV-2 virus. Similarly, [START_REF] Ito | Relative instantaneous reproduction number of Omicron SARS-CoV-2 variant with respect to the Delta variant in Denmark[END_REF] estimate that the effective reproduction number of the Omicron variant is more than 3 times that of the Delta variant in Denmark.

Our analysis of secondary case data did not, strictly speaking, allow one to conclude statistically that SSEs linked to the Delta variant had a different extreme value index from those linked to the original strains of SARS-CoV-2. However, in the contacttracing data recorded in South Korea, we did find a heavy tail in the offspring distribution when the Delta variant made the majority of cases, as opposed to when it did not. This tentative finding of a heavier tail in the data linked to the Delta

Figure 1 :

 1 Figure1: Top panels: Hill plots (solid lines) and corresponding 90% Gaussian asymptotic confidence intervals (dashed lines) as functions of the threshold value u, for n = 10,000 simulated data points Z i from the Burr distribution with probability density function f (x) = ξ -1 x -ρ/ξ-1 (1 + x -ρ/ξ ) 1/ρ-1 (for x > 0) with ξ = 1/2 and ρ = -1 in the left panel, and for the data Z i (i.e. the smallest integer larger than or equal to Z i ) in the right panel. Bottom panels: Averaged Hill plots when this experiment is repeated N = 1,000 times.

1 , 2 ,

 12 . . . with p ξ,σ (x) > 0 for ξ ≥ 0 or ξ < 0 and σ/ξ a negative integer, has been shown to outperform the continuous GPD when there are a large number of tied observations: see the simulated Poisson and discrete Inverse-Gamma examples in Section 3.1 of [13], which respectively show that the GPD provides poor model fits and poor tail estimates when the data is highly discrete, while the D-GPD distribution performs well. Its closed-form survival and probability mass functions allow for an exact likelihoodbased inference constructed upon the maximum likelihood estimators ξu , σu = arg min ξ>-1, σ>0 n i=1 log p ξ,σ (Z i -u).

  )) only partially contain the estimated extreme value index plot for SARS-CoV (solid blue line), reflecting a net difference between the two fat-tailed distributions of secondary cases associated with SARS-CoV and SARS-CoV-2. This conclusion is corroborated by the mean excess function estimates (Fig.2 (d)), which similarly indicate the relevance of separating the analysis for each coronavirus. This suggests that although SARS-CoV and SARS-CoV-2 belong to the same family of respiratory diseases, superspreading events are larger in scale for SARS-CoV in comparison to SARS-CoV-2. For all these reasons, pooling the data before applying extreme value tools can lead to misleading conclusions on the propagation of the SARS-CoV-2 virus.

  Negative binomial fits slope=1.43 => ξ=0.589 slope=0.256 => ξ=0.204 D-GPD log-mass functions

Figure 2 :

 2 Figure 2: Secondary case data from [3] (Datasets S1 and S2). (a) Histogram of the number of secondary cases for SARS-CoV (blue, n = 15) and SARS-CoV-2 (red, n = 45) SSEs. (b) Fitted probability mass function, conditional on Z > 6, of the negative binomial distribution for SARS-CoV (blue) and SARS-CoV-2 (red) SSEs.(c) Hill estimates of ξ for SSEs associated with SARS-CoV (solid blue), SARS-CoV-2 (solid red), and the pooled data (solid black), obtained from the exceedance values Z i -u given Z i ≥ u, as function of the threshold u, along with the resulting 90% confidence intervals for SARS-CoV (dashed blue) and SARS-CoV-2 (dashed red) SSEs. (d) Mean excess plots of SARS-CoV (blue) and SARS-CoV-2 (red) SSEs, quantified by the average of the exceedances Z i -u given Z i ≥ u, as function of u. (e) Discrete GPD maximum likelihood estimates of ξ for SARS-CoV (solid blue) and SARS-CoV-2 (solid red) SSEs, calculated from the exceedances Z i -u given Z i ≥ u, as function of u, along with their corresponding 90% confidence intervals (dashed lines), and the Hill plot produced by combining SARS-CoV and SARS-CoV-2 SSEs (black line). (f) Logarithm of the probability mass functions P σ,ξ (X = x) of the D-GPD fits to the exceedance values Z i -u given Z i ≥ u, for the thresholds u = 6 (dotted lines) and u = 10 (solid lines), for SARS-CoV (blue) and SARS-CoV-2 (red).

  3 (c)), and leads to an estimate of ξ around 0.24 with the 90% confidence intervals overwhelmingly suggesting an estimate greater than 0, thus confirming the fat-tailed nature of SARS-CoV-2 SSEs (Fig.3 (d)) in this sample. Interestingly, revisiting the small SARS-CoV-2 SSE dataset (Dataset S2) of size 45 using the D-GPD maximum likelihood estimation method (Fig.2(e)) results in an estimate of around 0.25, in agreement with the results from the Indian secondary case data. This suggests that the distribution of SARS-CoV-2 SSEs has a finite third moment and possibly even a fourth moment. These results are different from those obtained for the SARS-CoV SSEs. The latter rather point towards a distribution with infinite variance and thus a much heavier right tail. This is confirmed by considering the fitted D-GPD probability mass functions for secondary cases (Fig.2 (f)) that decrease much more rapidly for SARS-CoV-2 than for SARS-CoV.To examine the extreme value behaviour of the SARS-CoV-2 offspring distribution in different conditions, we turn to the analysis of two contact-tracing datasets in South Korea, a country which has a similar population density to the Indian state of Tamil Nadu, but did not resort to any full lockdown and has one of the largest and best-organised epidemic control programmes in the world. The first dataset was collected in the first half of 2020 (Database S4), while the second was collected during the fourth community epidemic in the summer of 2021 (Database S5) in the context of the assessment of transmission dynamics for the Delta variant of SARS-CoV-2.

  Tail index plot D-GPD -India data

Figure 3 :

 3 Figure 3: Secondary case data (Database S3) for SARS-CoV-2 from Andhra Pradesh and Tamil Nadu (India). (a) Bar plot of the log(Z i +1) (n = 88,527). (b) Mean excess plots of secondary cases. (c) Hill (solid black), continuous GPD maximum likelihood (solid blue) and discrete GPD maximum likelihood (solid bold red) estimates of ξ.(d) Discrete GPD maximum likelihood estimates of ξ (solid red) and their associated 90% confidence intervals (dashed red). In panels (c) and (d), the averaged discrete GPD estimate ξ = 0.239 over the stable region u ∈ [0, 10] is indicated with the horizontal red line.

Z

Fig. 5 (

 5 Fig.5(a) again backs the assumption of a heavy tail, but here, the mean excess plot in Fig.5(b) suggests a more convincing linearly increasing fit to the mean excess function with a slope of around 0.3. The Hill estimator and both continuous and discrete GPD maximum likelihood estimators, represented in Fig.5 (c), appear to support the fat tail assumption of the offspring distribution which is mainly dominated here by the Delta variant. Once again, the D-GPD estimate has a much smoother and more stable sample path, with a stable zone over u ∈[START_REF] Lloyd-Smith | Superspreading and the effect of individual variation on disease emergence[END_REF][START_REF] Swinkels | SARS-CoV-2 Superspreading Events Database[END_REF] indicating a point estimate of around 0.21. The 90% confidence interval of the D-GPD estimate over that region, provided in Fig.5 (d), does not contain 0 and offers further justification of the assumption that the offspring distribution is fat-tailed in this dataset, in contrast to the 2020 South Korea data where the validity of this conclusion is much less clear.

  Tail index plot D-GPD -Korea data

Figure 4 :

 4 Figure 4: Secondary case data (Database S4) for SARS-CoV-2 from South Korea (first half of 2020). (a) Bar plot of the log(Z i + 1) (n = 5,165). (b) Mean excess plots of secondary cases. (c) Hill (solid black), continuous GPD maximum likelihood (solid blue) and discrete GPD maximum likelihood (solid bold red) estimates of ξ. (d) Discrete GPD maximum likelihood estimates of ξ (solid red) and their associated 90% confidence intervals (dashed red).

  Tail index plot D-GPD -Delta data

Figure 5 :

 5 Figure 5: Secondary case data (Database S5) for SARS-CoV-2 from South Korea (July-August 2021). (a) Bar plot of the log(Z i + 1) (n = 33,903). (b) Mean excess plots of secondary cases. (c) Hill (solid black), continuous GPD maximum likelihood (solid blue) and discrete GPD maximum likelihood (solid bold red) estimates of ξ.(d) Discrete GPD maximum likelihood estimates of ξ (solid red) and their associated 90% confidence intervals (dashed red). In panels (c) and (d), the averaged discrete GPD estimate ξ = 0.209 over the stable region u ∈ [1, 10] is indicated with the horizontal red line.

Figs. 6

 6 Figs. 6 and 7 display the D-GPD maximum likelihood estimates of ξ as functions of the cluster size u. Eyeballed thresholds are indicated by the vertical dashed lines in Figs.6 and 7. The final selected estimates are reported in Table4, where[START_REF] Hitz | Discrete Extremes[END_REF]

6 : 7 :

 67 the estimated percentiles with actual observations to get a sense of what would constitute a conducive environment for the formation of large SARS-CoV-2 clusters.For example, the estimated 95th percentile of 120 cases in Kerala is close to two clusters of 113 cases (nursing home) and 132 cases (local transmission) already observed in Kerala. Likewise, the estimate q0.95 = 272 cases in Canada is fairly close to a cluster of 324 cases in Canadian nursing homes. In Oregon, the estimated 99th percentile q0.99 = 124 cases is in the vicinity of a cluster of 134 cases in a care home setting. In Colorado, the estimate q0.99 = 140 cases is close to a cluster of 134 cases in a nursing home. All of these clusters bar one (the local transmission cluster in Kerala) correspond to indoor environments where social distancing is difficult to practice. Analysis of cluster cases, for the four countries/states where the source is official data (Database S6). Plots of discrete GPD maximum likelihood estimates of ξ (solid lines), along with their 90% confidence intervals (dotted lines) and the final selected estimates (horizontal dashed lines) and thresholds (vertical dashed lines). Analysis of cluster cases as in Fig.

  during the early phase of the COVID-19 pandemic was fat-tailed, thus confirming their findings, although we show in various ways that this dataset should not be pooled with their 15 SSEs associated with SARS-CoV (Dataset S1), since they correspond to substantially different distributions. On the other hand, as accurate extreme value inference requires a large sample size in general, we also analysed an Indian secondary case dataset of size 88,527 collected in 2020 (Database S3), which contains a very large number of tied observations. The D-GPD estimate of the extreme value index is around 0.24, which is in full agreement with the estimate of around 0.25 found by revisiting the small SSE dataset of size 45 from[START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF]. The distribution of SARS-CoV-2 SSEs therefore appears to have at least a finite third moment, whereas that of SARS-CoV SSEs is found to have a much heavier upper tail with infinite variance and therefore stronger superspreading effect. In an effort to account for the quality of implemented control programmes as well as the nature of the variant under study, we used two extra South Korean contact-tracing datasets. For the first dataset (Database S4), collected in the first half of 2020 and used in[START_REF] Wong | Evidence that coronavirus superspreading is fat-tailed[END_REF], we cannot disprove that the distribution of the number of secondary cases has an exponential-type tail. By contrast, for the second South Korean dataset (Database S5) collected during the summer of 2021, in which the majority of cases correspond to the Delta variant of SARS-CoV-2[START_REF] Ryu | Serial interval and transmission dynamics during SARS-CoV-2 Delta variant predominance, South Korea[END_REF], we obtained a D-GPD estimate, ξ ≈ 0.21 clearly suggesting a heavier upper tail for the Delta variant and therefore more pronounced superspreading potential in South Korea relative to the first half of 2020.

  see for example the discussion in p.77 of Chapter 4 in[START_REF] Jacob | Extreme Value Statistics, in: Forecasting and Assessing Risk of Individual Electricity Peaks[END_REF]. It applies reasonably well to the D-GPD sample paths, because they are overall much smoother and more stable than the standard Hill and GPD maximum likelihood sample paths, which are not designed to handle the discreteness of the data. The development of more elaborate statistical techniques for the choice of threshold in discrete GPD maximum likelihood estimation, such as methods based on asymptotic MSE minimisation or the bootstrap in the spirit of the approaches outlined in Section 5.4 of[START_REF] Gomes | Extreme value theory and statistics of univariate extremes: A review[END_REF] for Hill estimation, is an open question which is beyond the scope of this paper.

Table 1 :

 1 Secondary case data (Database S3) for SARS-CoV-2 from Andhra Pradesh and Tamil Nadu (India).

	Z	0	1	2	3		4	5	6	7		8	9 10 11
	Count 62,540 17,493 4,885 1,730 802 444 267 149 67 44 29 22
	Z	12 13 14 15 16 17 18 19 21 22 23 25 28 31 37 39
	Count 14 16 3	3	4	4	1	1	2	1	1	1	1	1	1	1

Table 2 :

 2

		0	1		2	3	4 5 6 7 8 9
	Count 4,558 364 114 62 27 7 7 4 4 1
	Z	10 11 12 15 17 18 21 24 27 51
	Count 2	3	1	2	2	1	2	2	1	1

Secondary case data (Database S4) for SARS-CoV-2 collected in South Korea in the first half of 2020.

Table 3 :

 3 Secondary case data (Database S5) for SARS-CoV-2 collected in South Korea from 25th July 2021 to 15th August 2021.

	Z	0		1	2		3	4	5		6	7	8	9 10
	Count 29,193 2,154 1,121 594 332 207 113 53 53 21 21
	Z	11 12 13 14 15 16 17 18 19 21 22 24 32
	Count 6	8	5	3	3	2	2	3	3	1	2	2	1

Table 4 :

 4 Final results for SARS-CoV-2 cluster sizes by country (first column), the corresponding sample size n (second column), D-GPD maximum likelihood ξ estimate and 90% confidence interval (third column), selected cluster size threshold u and associated number n upon which the ξ estimate is calculated (fourth column), D-GPD maximum likelihood 95% and 99% percentile estimates of cluster size (fifth and sixth columns), and the sample maximum (last column). The top table corresponds to data from official sources (Database S6), and the bottom table to data from media sources (Database S7). The results reported in the latter table only concern the 9 countries and states for which the extreme value analysis was conclusive.

u of exceedance values Y i -u given Y i ≥ u

Funding. This research was supported by the French National Research Agency (grant numbers ANR-19-CE40-0013, ANR-17-EURE-0010), the TSE-HEC ACPR Chair and an AXA Research Fund Award on 'Mitigating risk in the wake of the COVID-19 pandemic'.

Acknowledgements. The authors acknowledge an anonymous Associate Editor and three anonymous reviewers for their helpful comments.

Australia

Competing interests. The authors declare no competing interests.