Simon Gabay

Ariane Pinche

Kelly Christensen

GALLIC(ORPOR)A: PROCESSING GALLICA'S HISTORICAL SOURCES

INTRODUCTION

The apparition of digital libraries have brought together unprecedented amounts of data and provided human sciences the means with which to integrate a new scientific paradigm. Such a transition implies overcoming a technological hurdle, particularly for historical documents which pose a significant number of challenges: manuscript documentation, difficult-to-read hands, variation-rich varieties of language. . . The Gallic(orpor)a project tries to overcome these obstacles. Because AI requires tremendous of data, the harmonisation of the production is a key problem.

TRAINING MODELS

HARMONISATION OF LAYOUT ANALYSIS

The SegmOnto controlled vocabulary is based on the assumption that most of the textual sources can be described in the same way if we use a codicological perspective focused on material aspects (running title, pagination/foliation, headings/rubrics, drop capitals, a main textual zone and potentially notes). This succeeds whether the sources are historical prints or manuscripts.

HARMONISATION OF TRANSCRIPTIONS

How to transcribe these two lines? ([. . . ]  apꝛes ꝯtraires toꝛmens  ml̕ t di / uers qil ot enu̕ s ymon lenchanteoꝛ). We offer recommendations regarding:

▶ Abbreviations such as ml̕ t for moult or enu̕ s for envers ▶ Tyronian notes such as  for et, or ꝯ for con.

▶ Allographetic (or graphic) variation such as vs s or ꝛ for r

DATA MODELLING

The information (text+layout) retrieved by HTR engines is encoded in ALTO-XML, which is not a format used by digital humanists. It needs to be converted in a standard format for text edition and publication: XML-TEI. Research had to be conducted on the optimal modelling, especially for the <sourceDoc> element.

<sourceDoc> <!--Page --> <surface xml:id="f1" n="0" cert="gold" ulx="0" uly="0" lrx="3800" lry="5600"> <graphic url="url"/> <!--TextBlock "MainZone:column#1"--> <zone xml:id="f1_z1" type="MainZone" subtype="column" n="1" points="800,580 800,4700 2100,4700 2100,580" source="url"> <!--TextLine "DefaultLine"--> <zone xml:id="f1_z1_l1" type="DefaultLine" subtype="subtype" n="1" points="1150,730 1140,640 2000,600 2030,700 2030,740 1150,759" source="url"> <!--baseline --> <path xml:id="f1_z1_l1_p" points="1150,730 2030,700"/> <!--transcription --> <line xml:id="f1_z1_l1_t">A transcription</line> </zone> </zone> </surface> </sourceDoc>

TOWARDS AUTOMATED CORPUS BUILDING

The advantage of the TEI is the extended flexibility and the richness of its <teiHeader> (i.e. metadata). Via API such as the BnF's SRU service, it is possible to automate the metadata retrieval along with the text contained in the digital facsimile to create automatically a corpus from a list of ARK identifiers (i.e. unique identifiers of the digital library).

FURTHER WORK

▶ Our models are not optimal for all the French sources: resources for incunabula, cursive writings, etc should be developed; ▶ Research is currently being made on the linguistic annotation of the data: lemmatisation, part-of-speech tagging, named entity recognition, linguistic normalisation. . . ; ▶ Producing data is not enough: an efficient web application should be offered to give access to the information; ▶ The philological exploitation of the data remains the main objective of this project. Books and manuscripts processed by our pipeline will be used to answer linguistic and linguistic questions using data-driven approach, such as stylometry or natural language processing.

  One of the main objective of the project is the creation of AI models to extract information from digital facsimiles. Because such a task is not trivial, empirical experiments are conducted to found the best configurations (network architectures, training data, scope of models. . . )

	Model	Accuracy
	Baseline	98.61%
	Gallic(orpor)a+	98.66%
	Gallic(orpor)a Antiqua	91.10%
	Gallic(orpor)a Antiqua/gothic	96.74%


Ian Hough
email: ian.hough@univ-grenoble-alpes.fr

Allan C Just

Bin Zhou

Michael Dorman

Johanna Lepeule

Itai Kloog

A multi-resolution air temperature model for France from MODIS and Landsat thermal data

Keywords: air temperature, land surface temperature, MODIS, Landsat, exposure error

Understanding and managing the health effects of ambient temperature (Ta) in a warming, urbanizing world requires spatially-and temporally-resolved Ta at high resolutions. This is challenging in a large area like France which includes highly variable topography, rural areas with few weather stations, and heterogeneous urban areas where Ta can vary at fine spatial scales. We have modeled daily Ta from 2000 -2016 at a base resolution of 1 km 2 across continental France and at a 200 x 200 m 2 resolution over large urban areas. For each day we predict three Ta measures: minimum (Tmin), mean (Tmean), and maximum (Tmax). We start by using linear mixed models to calibrate daily Ta observations from weather stations with remotely sensed MODIS land surface temperature (LST) and other spatial predictors (e.g. NDVI, elevation) on a 1 km 2 grid. We fill gaps where LST is missing (e.g. due to cloud cover) with additional mixed models that capture the relationship between predicted Ta at each location and observed Ta at nearby weather stations. The resulting 1 km Ta models perform very well, with ten-fold cross-validated R 2 of 0.92, 0.97, and 0.95, mean absolute error (MAE) of 1.4 °C, 0.9 °C, and 1.4 °C, and root mean square error (RMSE) of 1.9 °C, 1.3 °C, and 1.8 °C (Tmin, Tmean, and Tmax, respectively) for the initial calibration stage. To increase the spatial resolution over large urban areas, we train random forest and extreme gradient boosting models to predict the residuals (R) of the 1 km Ta predictions on a 200 x 200 m 2 grid. In this stage we replace MODIS LST and NDVI with composited top-of-atmosphere brightness temperature and NDVI from the Landsat 5, 7, and 8 satellites. We use a generalized additive model to ensemble the random forest and extreme gradient boosting predictions with weights that vary spatially and by the magnitude of the predicted residual. The 200 m models also perform well, with ten-fold cross-validated R 2 of 0. 79, 0.79, and 0.85, MAE of 0.4, 0.3, and 0.3, and RMSE of 0.6, 0.4, and 0.5 (Rmin, Rmean, and Rmax, respectively). Our model will reduce bias in epidemiological studies in France by improving Ta exposure assessment in both urban and rural areas, and our methodology demonstrates that MODIS and Landsat thermal data can be used to generate gapfree timeseries of daily minimum, maximum, and mean Ta at a 200 x 200 m 2 spatial resolution.

Introduction

Ambient or near-surface air temperature (Ta) is increasingly recognized as an important health risk. High or low Ta is associated with increased morbidity and mortality across regions and climates [START_REF] Gasparrini | Mortality risk attributable to high and low ambient temperature: A multicountry observational study[END_REF][START_REF] Guo | Global variation in the effects of ambient temperature on mortality: A systematic evaluation[END_REF][START_REF] Song | Impact of ambient temperature on morbidity and mortality: An overview of reviews[END_REF], and recent work suggests that high Ta may exacerbate the effect of exposure to particulate matter (PM), another major health hazard [START_REF] Li | Modification of the effects of air pollutants on mortality by temperature: A systematic review and meta-analysis[END_REF]. Ta exposure is a growing concern in cities, which are often warmer than the surrounding countryside due to increased heat accumulation and slower heat diffusion [START_REF] Arnfield | Two decades of urban climate research: A review of turbulence, exchanges of energy and water, and the urban heat island[END_REF]. Urban areas are now home to more than half the world's population, and this share is projected to increase to almost 70% by 2050 [START_REF] United Nations | World Urbanization Prospects: The 2018 Revision. USGS, 2018a. Landsat 4-7 Surface Reflectance (LEDAPS) Product Guide[END_REF]. Health effects of Ta are also seen in rural populations [START_REF] Lee | Study on the association between ambient temperature and mortality using spatially resolved exposure data[END_REF], although fewer studies have examined these due to the difficulty of estimating Ta exposure. Meanwhile climate change is increasing Ta and the frequency of extreme events such as heat waves in both urban and rural areas (IPCC, 2013). The health burden of Ta exposure is expected to grow as climate change and urbanization continue [START_REF] Gasparrini | Projections of temperature-related excess mortality under climate change scenarios[END_REF][START_REF] Wang | Accounting for adaptation and intensity in projecting heat wave-related mortality[END_REF].

Understanding, monitoring, and managing Ta health effects requires spatiotemporally-resolved Ta at high resolutions. Weather station networks measure Ta at high temporal resolution, but rarely capture spatial variation at the scales needed for epidemiological studies (e.g. across a region, within a city). Failure to account for spatial variation in Ta can introduce error in exposure assessment, which tends to bias health effect estimates towards the null [START_REF] Zeger | Exposure measurement error in time-series studies of air pollution: Concepts and consequences[END_REF]. Some recent epidemiological studies have addressed this issue by using spatiotemporally-resolved Ta estimates from numerical weather prediction models such as WRF (Ha et al., 2017b(Ha et al., , 2017a)), but computational limitations currently restrict these models to medium spatial resolutions (e.g. 4 km) or small geographic areas (e.g. a single city). In urban areas, studies have used weather model Ta estimates or indicators such as sky view factor, vegetation abundance, and land surface temperature to create indexes that identify warmer and cooler areas within a city [START_REF] Goggins | Effect modification of the association between short-term meteorological factors and mortality by urban heat islands in Hong Kong[END_REF][START_REF] Ho | Delineation of spatial variability in the temperature-mortality relationship on extremely hot days in greater Vancouver, Canada[END_REF][START_REF] Laaidi | The impact of heat islands on mortality in Paris during the August 2003 heat wave[END_REF][START_REF] Milojevic | Methods to estimate acclimatization to urban heat island effects on heat-and cold-related mortality[END_REF][START_REF] Smargiassi | Variation of daily warm season mortality as a function of micro-urban heat islands[END_REF]. Studies to date have focused on the typical spatial distribution of Ta during a specific time period (e.g. a single heat wave, the hot season) as the limited temporal variability of the indicator variables and cost of numerical weather prediction have precluded consideration of changes in the pattern of warmer and cooler areas over time.

Other recent studies have used Ta estimates from hybrid land use regression models that predict Ta based on remotely sensed 1 km land surface temperature (LST) and spatial and spatiotemporal variables such as elevation and normalized difference vegetation index (NDVI) [START_REF] Kloog | Using Satellite-Based Spatiotemporal Resolved Air Temperature Exposure to Study the Association between Ambient Air Temperature and Birth Outcomes in Massachusetts[END_REF]Shi et al., 2016b[START_REF] Shi | Impacts of temperature and its variability on mortality in New England[END_REF]. This approach takes advantage of the growing body of satellite earth observation data and the fact that LST is a good indicator of spatiotemporal variation in Ta [START_REF] Oyler | Remotely sensed land skin temperature as a spatial predictor of air temperature across the conterminous United States[END_REF]. In particular, a technique that uses linear mixed models to calibrate the relationship between daily 1 km LST from the Moderate Resolution Imaging Spectroradiometer (MODIS) instrument and Ta has been shown to perform well over large, heterogeneous areas including the northeastern USA (root mean square error

[RMSE] 2.2 °C) [START_REF] Kloog | Predicting spatiotemporal mean air temperature using MODIS satellite surface temperature measurements across the Northeastern USA[END_REF], the southeastern USA (RMSE 1.4 °C) (Shi et al., 2016a), France (RMSE 1.7 °C) [START_REF] Kloog | Modelling spatio-temporally resolved air temperature across the complex geo-climate area of France using satellite-derived land surface temperature data[END_REF], and Israel (RMSE 1.2 °C) [START_REF] Rosenfeld | Estimating daily minimum, maximum, and mean near surface air temperature using hybrid satellite models across Israel[END_REF].

These models are parsimonious compared to numerical weather prediction, which allows them to capture both spatial and temporal variation in Ta over large areas and long time periods. Their spatial resolution suffices for areas where Ta varies little at scales of less than 1 km and for studies where subjects' locations are only approximately known. But finer spatial resolution estimates are needed for studies with address-level location data, particularly in urban areas where Ta can vary markedly within a square kilometer. Very high spatiotemporal resolutions would also benefit studies that have time-location data (e.g. GPS tracks).

In this study we extend the mixed modeling approach to predict daily minimum, mean, and maximum Ta (Tmin, Tmean, Tmax, respectively) at a 1 km resolution across continental France and at a 200 m resolution across 103 urban areas in continental France. We improve performance at the 1 km resolution by allowing the daily Ta ~ LST relationship to vary between climatic regions, and we consider both daytime and nighttime MODIS LST, which allows us to predict diurnal (Tmax) and nocturnal (Tmin) temperature in addition to Tmean. This is useful both for studies of urban heat islands, which exhibit different spatial patterns and intensities during day vs. night [START_REF] Arnfield | Two decades of urban climate research: A review of turbulence, exchanges of energy and water, and the urban heat island[END_REF], and for studies of Ta variability, which recent work suggests may independently affect health [START_REF] Guo | Temperature variability and mortality: A multi-country study[END_REF][START_REF] Molina | The perils of climate change: In utero exposure to temperature variability and birth outcomes in the Andean region[END_REF][START_REF] Shi | Impacts of temperature and its variability on mortality in New England[END_REF].

We also add a local stage that uses an ensemble of machine learning algorithms to predict the residuals of the 1 km model in urban areas based on higher spatial resolution predictors including thermal data from the Landsat 5, 7, and 8 satellites. This allows us to predict daily Ta over 17 years at a 200 m spatial resolution which better captures intra-urban Ta variation across 103 urban areas.

Data and methods

Study area and period

Our study area is continental France, comprising all French territory in Europe except Corsica.

It covers 542,973 km 2 of topographically and climatically diverse terrain with elevations that range from -10 to 4,809 m. [START_REF] Joly | Les types de climats en France, une construction spatiale[END_REF] classify France into eight climatic regions based on the magnitude, variability, and seasonality of temperature and precipitation (Fig. 1). The north and west coasts have a wet, temperate oceanic climate, which transitions to a drier, cooler modified oceanic climate in the north center. The mountainous east, south center, and southwest have variable montane and semi-continental climates with cold winters. In the southeast, the Mediterranean coast has hot, dry summers with mild wet winters; the inland southeast and isolated segments of the west coast are similar but cooler. The southwest basin resembles the inland southeast but with drier winters. The estimated population on January 1, 2018was 64,388,583 (INSEE, 2018). About 80% of the population is urban, and this share is projected to grow to 88% by 2050 [START_REF] United Nations | World Urbanization Prospects: The 2018 Revision. USGS, 2018a. Landsat 4-7 Surface Reflectance (LEDAPS) Product Guide[END_REF]. The largest urban area, Paris, has a population of 12.5 million (20% of the total) and the six next largest urban areas have a population of one to 2.3 million (combined 14% of total). A further 10% of the population lives in cities of one half to one million, and 37% live in urban areas with fewer than half a million residents (Fig. S1). Our study period is January 1, 2000 through December 31, 2016.

Meteorological observations

We use daily weather station observations from Météo France, the French national meteorological service. About 64% of the observations come from stations managed by Météo France; the remaining stations are managed by other entities. All observations are quality controlled by Météo France. We exclude stations with no metadata or that do not record hourly Ta, and for each month during the study period we exclude stations that were active for fewer than 21 days in the month. This leaves 1710 to 2314 stations on each day. The stations are distributed over the entire study region, but are denser in populous areas (e.g. Paris, the southeast) and the Alps (which has many ski resorts, hydroelectric dams, and avalanche monitors) (Fig. 1). Just 3% of the stations are located within large urban areas (as defined in section 2.7), 7% are in peri-urban areas (within 5 km of an urban area), and the remaining 90% are rural.

The stations calculate daily Tmin as the lowest Ta observed from 18 UTC the previous day until 18 UTC on the day; daily Tmax is the highest Ta observed from 6:00 UTC on the day until 6:00 UTC the following day. Most stations calculate Tmean as the mean of all (at least 24) Ta observations from 0 UTC on the day until 0 UTC the following day. However, about 40% of the Tmean observations were calculated as the average of Tmin and Tmax. We exclude these observations, meaning our final dataset has fewer observations for Tmean than for Tmin or Tmax.

Daily Ta at the included stations during the study period ranged from Tmin of -31.2 °C to Tmax of 44.1 °C; mean Tmean was 11.3 °C with a standard deviation of 7.1 °C (Table S1).

Land surface temperature and emissivity

We use version 6 of the widely-used MODIS daily 1 km land surface temperature and emissivity product from the Terra and Aqua satellites (MOD11A1 and MYD11A1, respectively) (Table 1). These products include daytime and nighttime LST derived using a split-window algorithm and land use classification-based emissivity and have been masked for clouds and validated to ± 2 K in clear-sky conditions across 47 sites on all seven continents [START_REF] Wan | New refinements and validation of the collection-6 MODIS land-surface temperature/emissivity product[END_REF]. We use the quality assessment band to exclude pixels with LST error > 2 K. As LST retrieval error increases over snow and water, we also exclude pixels with NDVI < 0 or where the corresponding 1 km grid cell has land cover of > 33% water.

Top-of-atmosphere brightness temperature

For large urban areas, we composite 30 m top-of-atmosphere brightness temperature (Tb) from the Landsat 5, 7, and 8 satellites (Table 1). Tb is the kinetic temperature a perfect blackbody would have if it emitted the quantity of thermal radiation measured by the satellite instrument.

Converting Tb to LST requires correcting for atmospheric effects and accounting for the emissivity of the earth's surface. This is difficult in the case of the Landsat satellites because Landsat 5 and 7 have only a single thermal band and the USGS Landsat 8's second thermal band is contaminated by stray light, precluding the use of the split-window algorithm [START_REF] Li | Satellitederived land surface temperature: Current status and perspectives[END_REF]. A global Landsat LST product is under development but not yet available [START_REF] Malakar | An Operational Land Surface Temperature Product for Landsat Thermal Data: Methodology and Validation[END_REF], so for this study we use Tb from the USGS Landsat Collection 1 Level-2 surface reflectance products (USGS, 2018a(USGS, , 2018b)).

The 16-day revisit time of the Landsat satellites means that Tb is unavailable for many locations on many days. Cloud cover and sensor malfunctions also contribute to these data gaps and can increase error in Tb retrieval. To reduce error, we discard all scenes with cloud cover > 75%.

We also discard all scenes captured during periods of instrument malfunction, which we identified by checking summary statistics of each scene for unrealistic values (e.g. mean Tb > 100 °C). We then trim the edges of Landsat 5 scenes by 2.5 km to remove abnormalities [START_REF] Robinson | A dynamic landsat derived normalized difference vegetation index (NDVI) product for the conterminous United States[END_REF] and mask pixels identified as high-or medium-confidence cloud in the pixel quality assessment band. We mask any remaining pixels where Tb ≤ -25 °C or Tb ≥ 50 °C.

Finally, for each calendar month we composite all Tb retrievals during the entire study period (e.g. every January in 2000 -2016). This yields 12 gap-free Tb datasets representing the 17year mean Tb of each pixel in each calendar month.

NDVI

We use version 6 of the MODIS monthly composite 1 km NDVI product from the Terra and Aqua satellites (MOD13A3 and MYD13A3, respectively). For large urban areas we also composite 30 m NDVI from the Landsat 5, 7, and 8 Collection 1 Level-2 surface reflectance products. We use a similar quality assurance and compositing procedure as for Tb, first discarding all scenes with greater than 75% cloud cover or during periods of thermal sensor malfunction (as this results in unreliable cloud confidence scores in the pixel quality assessment band). We then trim the edges of Landsat 5 scenes by 2.5 km and adjust NDVI from Landsat 5 and Landsat 7 to match Landsat 8 using equation Eq. 1 [START_REF] Robinson | A dynamic landsat derived normalized difference vegetation index (NDVI) product for the conterminous United States[END_REF].

NDVI !" = 0.0235 + 0.9723 × NDVI #$,#&
Eq. 1

Similar to [START_REF] Robinson | A dynamic landsat derived normalized difference vegetation index (NDVI) product for the conterminous United States[END_REF], for each calendar month we create two 17-year mean composites, one using pixels marked as clear in the pixel quality assurance band (i.e. not cloud, cloud shadow, snow, or water) and a second using pixels marked as snow or water. Finally, we mosaic the two composites preferring the clear pixels composite.

Elevation, Population, Land Cover, and Climatic Regions

We thematic accuracy [START_REF] Bossard | CORINE land cover technical guide -Addendum[END_REF]. We aggregate the land cover classes into four groups:

artificial, vegetation, bare, and water (Table S2). Finally, we use the eight climatic regions of [START_REF] Joly | Les types de climats en France, une construction spatiale[END_REF], which are based on temperature and precipitation patterns (Fig. 1).

Model grids

For the 1 km model, we create a grid covering continental France by making a 1 km square buffer around the centroid of each MODIS 1 km LST pixel in the ETRS89-LAEA Europe (EPSG:3035) equal-area projection. We associate each 1 km grid cell with the MODIS LST and NDVI pixel having the same centroid and calculate the mean elevation, total population, percent area of each land cover group, and climate region with greatest spatial overlap.

For the 200 m model, we create a grid covering large urban areas. Starting from a 200 m grid in the ETRS89-LAEA Europe (EPSG:3035) equal-area projection, we select all cells in continental France containing "Urban fabric" or "Industrial or commercial units" in the 2012 CLC inventory. We associate each cell with the corresponding INSEE gridded population and select cells with 50 or more inhabitants as well as the eight surrounding cells (i.e. including diagonal neighbors). We define urban areas as four-wise contiguous (i.e. excluding diagonal neighbors) groups of cells and sum the population of all cells in each urban area. Finally, we eliminate urban areas with population < 50,000. This leaves 103 large urban areas ranging from greater Paris (9.4 million inhabitants) to Armentières (50,260 inhabitants). For each 200 m grid cell in a large urban area or that contains a weather station we calculate the mean 17-year composite Landsat Tb and NDVI for each calendar month, mean elevation, and percent area of each land cover group.

Statistical methods

We use a four-stage approach to predict Ta: stages 1 and 2 predict daily 1 km Ta across continental France and stages 3 and 4 predict daily 200 m Ta within large urban areas. We consider each year during the study period (2000 -2016) and each Ta measure (Tmin, Tmax, and Tmean) separately. Stages 1 and 2 are an extension of the method used in [START_REF] Kloog | Modelling spatio-temporally resolved air temperature across the complex geo-climate area of France using satellite-derived land surface temperature data[END_REF] and are detailed in Appendix A. Sections 2.8.1 to 2.8.2 detail stages 3 and 4; the following is a brief overview of all stages.

In stage 1 we calibrate Ta at each station as a function of daily 1 km LST and emissivity, monthly 1 km NDVI, and 1 km elevation, population, and land cover. We use a linear mixed model to allow the Ta ~ LST relationship to vary by day within each climatic region. We use this calibrated relationship to predict 1 km Ta (Tap_s1) for all cell-days where LST is available.

In stage 2, we fill gaps in Tap_s1 where 1 km LST was not available by calibrating Tap_s1 as a function of daily 1 km inverse distance weighting interpolated observed Ta (TIDW). We use a linear mixed model to allow the Tap_s1 ~ TIDW relationship to vary by location. We use this calibrated relationship to fill gaps in Tap_s1, producing gap-free daily 1 km predicted Ta (Tap_1km). This is the 1 km Ta model.

In In stage 4, we calibrate a generalized additive model that ensembles Rp_rf and Rp_gb. We use a tensor product smooth with interaction to allow the relative performance of the RF and GB models to vary by location and with the magnitude of the predicted residual. Finally, we add the ensemble predictions to Tap_1km to get daily 200 m predicted Ta for large urban areas (Tap_200m). This is the 200 m Ta model. Populationi is the population of cell i; xi and yi are the geographical coordinates of cell i; j is the Julian day; and eij is the error for cell i on day j.

We use the R packages ranger [START_REF] Wright | Ranger: A fast implementation of random forests for high dimensional data in C++ and R[END_REF], XGBoost (Chen and Guestrin, 2016), and mlr [START_REF] Bischl | Mlr: Machine learning in R[END_REF] to train the random forest and XGBoost models. We tune the models using the sequential model-based optimization of package mlrMBO [START_REF] Bischl | mlrMBO: A Modular Framework for Model-Based Optimization of Expensive Black-Box Functions[END_REF]. Briefly, mlrMBO estimates optimal hyperparameter values by iteratively training and evaluating a model using hyperparameter values that are chosen based on the performance of previous iterations. We use a fixed number of iterations and evaluate performance as the mean RMSE of two random 80% holdouts (i.e. we train the model on a 20% random sample of the data, predict and calculate RMSE for the held-out 80%, repeat, and take the mean of the two RMSEs). Initial exploration showed that this resampling approach produced stable estimates of RMSE at a lower computational cost than cross-validation.

For the random forest, we use 400 trees and a minimum of 5 observations per node, and tune mtry (the number of variables to consider for each split) from 3 to 12 (25% to 100% of the explanatory variables) using 6 mlrMBO iterations. Initial exploration showed that using more than 400 trees only marginally increased performance and had a high computational cost. For the XGBoost model, we use the gbtree booster with 100 rounds and set gamma (the minimum loss reduction for a split) to 5. We use 24 mlrMBO iterations to tune eta (the learning rate) from 0.1 to 0.3, the maximum tree depth from 5 to 20, the minimum number of observations per node from 3 to 30, and the fraction of features used in each tree from 0.75 to 1.

We evaluate the performance of the stage 3 models using 5-fold cross-validation with nested tuning. We use the final stage 3 random forest and XGBoost models to predict the residual of the 1 km Ta model (Rp_rf and Rp_xgb, respectively) for all 200 m cell-days.

Performance assessment

We use 10-fold out-of-sample cross-validation to assess the overall performance of the models.

For the random forest and XGBoost model we use nested tuning (i.e. within each crossvalidation fold we tune the model as described in section 2.8.1). To evaluate the models' ability to capture both spatial and temporal patterns in Ta, we also calculate the spatial and temporal components of the errors. The spatial component is the difference at each station between the annual mean of daily observed Ta (T ) ), and the annual mean of daily predicted Ta (T )* ). The temporal component is the difference at each station between DTa and DTap where DTa is the difference between daily observed Ta and T ) and DTap is the difference between daily predicted Ta and T )* . We use Google Earth Engine [START_REF] Gorelick | Google Earth Engine: Planetary-scale geospatial analysis for everyone[END_REF] to quality assure and composite Landsat Tb and NDVI and aggregate them to the 200 m grid cells. For all other data processing and analyses we use R version 3.4.4 (R Core Team, 2018).

Results

Table 2 presents the mean 10-fold cross-validated performance of the stage 1 models (predicting daily 1 km Ta from LST) across all years. The models perform very well, with R 2 of 0.92 or higher, RMSE of less than 2 °C, and mean absolute error (MAE) of less than 1.5 °C.

All models have very low bias: the slope of observed vs. predicted Ta is 1.00 while the intercept ranges from 0.01 to 0.02. The Tmean models perform best overall (MAE 0.94), followed by the Tmax (MAE 1.35) and Tmin (MAE 1.43) models. The models capture both spatial and temporal variation in Ta and show little variation in performance between years, although overall Tmean performance decreases slightly after 2010, possibly reflecting degradation of the Terra MODIS instrument (Table S4). Consistent with previous studies, nighttime LST is the best predictor of Tmin and Tmean while daytime LST is the best predictor of Tmax [START_REF] Oyler | Remotely sensed land skin temperature as a spatial predictor of air temperature across the conterminous United States[END_REF][START_REF] Rosenfeld | Estimating daily minimum, maximum, and mean near surface air temperature using hybrid satellite models across Israel[END_REF][START_REF] Yoo | Estimation of daily maximum and minimum air temperatures in urban landscapes using MODIS time series satellite data[END_REF]. Aqua LST is a better predictor of Tmin and Tmax while Terra LST is a better predictor of Tmean. This is expected as the Aqua overpasses (approximately 1:30 and 13:30 local solar time) are closer to the time at which Tmin and Tmax typically occur in France.

However, Aqua LST is only available since July 2002, so we use Terra LST for all models prior to 2003. 4 presents the performance by climatic region and urban vs. rural locations. The Tmin and Tmean models perform slightly less well in winter months, possibly due to higher LST missingness from more frequent cloud cover. The Tmax model performs best in late winter, early spring, and fall. The models perform less well in the mountain, semi-continental, and modified Mediterranean climates. These climates occur in mountainous areas where large contrasts in topography and land cover make modelling particularly challenging; other factors not included in the model may also reduce performance in these areas. The models perform slightly better in peri-urban areas than in urban and rural areas, possibly due to the higher density of weather stations (peri-urban areas have the most stations per km 2 ). On the hot summer day of Aug 10, 2012, predictions ranged from a Tmin of 3 °C in parts of the Alps to a Tmax of 39 °C in the southeast and southwest. On the Tmin map, the southwestern cities of Toulouse and Bordeaux stand out as hotspots, while Paris and Rouen are faintly visible as warm spots in the north. The north is colder than the Vosges mountains in the northeast and the Pyrenees in the southwest are warmer than the alps. The warmest areas are the southern Rhone river valley in the southeast and a patch of the southwestern Atlantic coast. On the Tmean map, Paris and Rouen are still visible, Lyon stands out in the east, and a few northwestern cities appear. Much of the southwest is as warm as the southeast, and the southwestern cities are harder to distinguish from the countryside. On the Tmax map, Lyon, Rouen, and some northwestern cities remain faintly visible, Pau and Tarbes appear in the southwest, and the north is warmer than the Vosges. 6 presents the performance by climatic region and urban vs. rural locations. These models also perform well, with overall R 2 of 0.79 to 0.85, RMSE of 0.41 to 0.63, and MAE of 0.26 to 0.39 (residual scale). As with the stage 1 models, the RTmean predictions are slightly better than the RTmin or RTmax predictions and the models perform least well in the mountain, semi-continental, and modified Mediterranean climates. The RTmin model performs slightly worse in late summer; otherwise performance is quite consistent across months and seasons.

The models have low bias, with a slope of observed vs. predicted of 1.00 and intercept of zero for every year. Performance is consistent across years except for the RTmin model, which performs slightly better in 2000 -2002, and the RTmean model, which performs best in 2004 (Table S6). Spatial location and elevation are generally the most important features in the random forest and XGBoost models (Fig. S2 -S3). Day of year and predicted 1 km Ta were equally or even more important in some models but less important in others. Landsat Tb and NDVI and population also contributed to the models, particularly for RTmean. The land cover and climatic region variables were the least important. Tmin was cool and capture the Garonne river in the center. The northwestern and northeastern suburbs have greater contrast with some areas slightly cooler than in the 1 km predictions and others slightly warmer. In the small city of Nancy, at 1 km both the city center and an area of ponds to the southeast have Tmin about 2 °C warmer than the surroundings. The 200 m predictions show warmer Tmin throughout most of the built-up area with sharp contrasts between built and open areas: compared to the 1 km predictions, Tmin is up to 2 °C higher in the center, north, and west of the built-up area and up to 2 °C lower over parks and over fields abutting the eastern edges of the city.

Discussion

Spatiotemporally-resolved Ta at high resolutions is essential to understanding, monitoring, and managing the health effects of Ta, a pressing issue in a warming, urbanizing world. We have developed the longest (2000 -2016), highest spatial resolution (1 km) model of daily Ta available for continental France aimed at public health research. Furthermore, our model provides an unprecedented spatial resolution of 200 m over large urban areas.

A key feature of our model is its ability to capture spatial variation in Ta. Previous epidemiological research in France linked geographical variation in mortality risk to both typical [START_REF] Laaidi | Temperature-related mortality in France, a comparison between regions with different climates from the perspective of global warming[END_REF] and extreme Ta [START_REF] Le Tertre | Impact of the 2003 heatwave on all-cause mortality in 9 French cities[END_REF]) using weather stations.

Recent studies in the USA showed that a daily 1 km Ta dataset similar to ours was needed to detect associations with low birth weight [START_REF] Kloog | Using Satellite-Based Spatiotemporal Resolved Air Temperature Exposure to Study the Association between Ambient Air Temperature and Birth Outcomes in Massachusetts[END_REF] and mortality [START_REF] Shi | Impacts of temperature and its variability on mortality in New England[END_REF]. Our model will allow future studies in France to include participants in rural areas far from weather stations and will also improve exposure estimates in urban areas.

Another key feature is our model's 200 m spatial resolution over urban areas. Estimating Ta exposure in cities is particularly challenging due to complex built environments and the scarcity of representative Ta measurements, as weather stations tend to be located outside cities (e.g. at airports) or in large parks. Consequently, few epidemiological studies have examined intraurban variation in Ta. In Milan, Italy, de'Donato et al. (2008) found that on hot summer days temperature measured at a nearby airport tended to be higher and more strongly associated with mortality than temperature measured in the city center, but in Turin and Rome there was little difference in temperature or its association with mortality between the city center and a nearby airport. In Paris, France, [START_REF] Laaidi | The impact of heat islands on mortality in Paris during the August 2003 heat wave[END_REF] used 1 km LST as a proxy for Ta and found an association between minimum LST and mortality during the August 2003 heatwave. In Brisbane, Australia, [START_REF] Guo | Spatiotemporal model or time series model for assessing citywide temperature effects on mortality?[END_REF] found no significant difference in the mortality ~ Ta relationship when estimating Ta exposure using a central weather station vs. kriging, but noted that there was little spatial variation in temperature across the city. In Seattle, USA, [START_REF] Ho | Delineation of spatial variability in the temperature-mortality relationship on extremely hot days in greater Vancouver, Canada[END_REF] found a significant association between spatial variation in mortality on extremely hot days and modeled humidex (a measure of both Ta and humidity). Our model will help future studies clarify the health effects of intra-urban Ta variation. Our model's unique combination of lower spatial resolution (1 km) predictions over a large geographical extent and higher spatial resolution (200 m) predictions over more densely populated areas will be particularly helpful for epidemiological studies. Broad geographical coverage is essential to including rural residents which have often been excluded from epidemiological studies, especially in France where the 103 largest urban areas covered by our 200 m Ta model contain less than half of the population. At the same time, high spatial resolution is important in dense urban areas where Ta can vary at fine spatial scales and the effect of spatial Ta variation is less well understood. Limiting the 200 m resolution predictions to large urban areas reduces computational effort while still covering a large portion of the population.

A fourth feature of our model is its ability to predict daily Tmin, Tmean, and Tmax. While Tmean suffices for many health studies [START_REF] Barnett | What measure of temperature is the best predictor of mortality?[END_REF], certain research questions may benefit from having Tmin and Tmax. For example, heatwave studies may wish to use heatwave definitions that refer to Tmin or Tmax [START_REF] Xu | Impact of heatwave on mortality under different heatwave definitions: A systematic review and meta-analysis[END_REF] or explore whether certain populations are sensitive to Tmin or nighttime Ta [START_REF] Laaidi | The impact of heat islands on mortality in Paris during the August 2003 heat wave[END_REF][START_REF] Murage | Effect of night-time temperatures on cause and age-specific mortality in London[END_REF]. Tmax might also be of interest because it tends to occur in the afternoon when people are more likely to be outside and active [START_REF] Guo | Heat wave and mortality: A multicountry, multicommunity study[END_REF]. Tmin and Tmax also allow calculating diurnal Ta range for studies of Ta variability and delineating diurnal and nocturnal urban heat islands for urban climate studies.

We demonstrate that allowing the relationship between 1 km LST and Ta to vary by climatic region as well as by day slightly improves performance: our stage 1 Tmean model achieves overall R 2 of 0.97 with RMSE of 1.29 whereas an initial version achieved R 2 of 0.96 with RMSE of 1.52 [START_REF] Kloog | Modelling spatio-temporally resolved air temperature across the complex geo-climate area of France using satellite-derived land surface temperature data[END_REF]. We also demonstrate that a GAM ensemble of machine learning models can use higher spatial resolution predictors including Landsat thermal data to account for some of the residual error in our daily 1 km Ta predictions. Adding this local stage both increases the spatial resolution of our model and improves performance.

One limitation of our method is its reliance on historical satellite thermal data. Our model is restricted to the MODIS period of record, which starts in 2000. Older thermal data is available from other satellites (e.g. Landsat), but not with a twice-daily revisit time. In the USA, [START_REF] Oyler | Creating a topoclimatic daily air temperature dataset for the conterminous United States using homogenized station data and remotely sensed land skin temperature[END_REF] showed that an anomaly-climatology approach could model daily Tmin and Tmax since 1948 from 8-day composite MODIS LST, although their approach may have smoothed spatiotemporal Ta trends.

Our model can estimate past Ta but, unlike numerical weather prediction models, cannot forecast future Ta. However, our model is much simpler, which allows us to run it at relatively high spatial resolutions (1 km and 200 m). In comparison, Météo France's weather prediction model has run at a spatial resolution of 1.3 km only since 2015, and the ECMWF's most recent ERA5 reanalysis has a spatial resolution of just 30 km. And recent studies suggest that incorporating LST from geostationary satellites might allow us to estimate close to real-time Ta [START_REF] Bechtel | Towards a satellite based monitoring of urban air temperatures[END_REF][START_REF] Keramitsoglou | An online system for nowcasting satellite derived temperatures for urban areas[END_REF], or possibly forecast next-day Ta from presentday MODIS LST [START_REF] Yoo | Estimation of daily maximum and minimum air temperatures in urban landscapes using MODIS time series satellite data[END_REF].

Another limitation of our approach is the temporal misalignment between observations of LST and Ta in the stage 1 model: the satellite overpass does not always coincide with the time that Tmin or Tmax occurs. Our model's low MAE (typically less than 1.5 °C) suggests that it produces good Ta estimates despite this; incorporating high temporal-resolution (e.g. hourly) LST from geostationary satellites might improve performance.

A fourth limitation of our model is the need to fill gaps in satellite thermal data. This can introduce error and may make modelling impossible in some areas or time periods. Landsat data is particularly challenging due to the satellites' 16-day revisit time; parts of France have no usable Landsat observations during some winters. The few previous studies that used Landsat thermal data to model Ta limited their analysis to days and locations where Landsat data was available [START_REF] Pelta | Spatiotemporal estimation of air temperature patterns at the street level using high resolution satellite imagery[END_REF] or used a few scenes that were deemed typical of hot summer days [START_REF] Ho | A comparison of urban heat islands mapped using skin temperature, air temperature, and apparent temperature (Humidex), for the greater Vancouver area[END_REF][START_REF] Ho | Mapping maximum urban air temperature on hot summer days[END_REF][START_REF] Wicki | Evaluation and Modeling of Urban Heat Island Intensity in Basel, Switzerland[END_REF]. We fill gaps in Landsat Tb by compositing all scenes for each calendar month across 17 years. This smooths spatial patterns and means we rely entirely on MODIS to capture short-term temporal variation in LST.

Combining data from Landsat 5, 7, and 8 may also introduce error as the sensors operate at different wavelengths and spatial resolutions (Table 1). Future studies may benefit from the forthcoming Landsat Surface Temperature product [START_REF] Malakar | An Operational Land Surface Temperature Product for Landsat Thermal Data: Methodology and Validation[END_REF] which might be more consistent, and would allow using LST as a predictor rather than brightness temperature.

Future studies could also make use of high spatial-resolution LST from forthcoming satellites. Landsat 9 will have a spatial resolution and revisit time similar to the previous Landsat satellites, but should offer better LST retrieval thanks to the correction of the stray light contamination that affects Landsat 8 [START_REF] Hair | Landsat 9 thermal infrared sensor 2 architecture and design[END_REF]. HyspIRI aims to provide a 60 m spatial resolution with a revisit time of 5 days [START_REF] Lee | An introduction to the NASA Hyperspectral InfraRed Imager (HyspIRI) mission and preparatory activities[END_REF], while MISTIGRI aims for 50 m spatial resolution with a daily revisit, but with coverage only within 15 ground tracks [START_REF] Lagouarde | The MISTIGRI thermal infrared project: Scientific objectives and mission specifications[END_REF]. If these satellites improve LST retrieval and reduce missingness then they could improve our method's ability to capture Ta over urban areas.

MODIS LST also contains gaps, which we do not fill. Rather, we predict daily 1 km Ta only where MODIS LST is available and fill gaps in the predictions based on nearby Ta observations. [START_REF] Li | Developing a 1 km resolution daily air temperature dataset for urban and surrounding areas in the conterminous United States[END_REF] achieved similar performance (RMSE 2.1 °C Tmin, 1.9 °C Tmax) for urban and surrounding areas in the USA by first filling gaps in MODIS LST using spatiotemporally nearby LST observations and then predicting daily Ta using geographically weighted regression. These approaches both assume that the spatial distribution of Ta or LST is similar on clear and cloudy days. [START_REF] Zhu | Retrievals of all-weather daytime air temperature from MODIS products[END_REF] used the MODIS atmospheric profile and cloud cover products to estimate instantaneous Ta in parts of China and the USA. Their approach had the additional advantage of not requiring any weather station Ta observations to calibrate the model, but it produced larger errors (RMSE 3.4 °C China, 2.9, USA). Despite these limitations, our model provides very good predictions of historical daily Ta for continental France at a 1 km or finer spatial resolution. These predictions may help compare rural and urban populations, identify and monitor urban heat islands, and better understand health effects. More broadly, our methodology and predictions may be useful in other geographical areas and for any application where Ta is a key variable.

Declarations of interest

None.

Fig. 1 .

 1 Fig. 1. Climatic regions of France according to Joly et al. (2010) and METEO-FRANCE stations used in the current study.

  use version 1.1 of the European Digital Elevation Model (EU-DEM) from the Copernicus Land Monitoring Service. These data have a 25 m spatial resolution and vertical RMSE of ±7 m (Tøttrup, 2014). We also use 200 m gridded 2010 population from INSEE, the French national statistics agency (INSEE, 2017). We use the 100 m Corine Land Cover (CLC) inventory for 2000, 2006, and 2012. The 2000 edition has been validated to better than 85%

  stage 3, we calculate the daily 200 m residuals of the 1 km Ta model (R) and train random forest (RF) and extreme gradient boosting (GB) models to predict R based on latitude, longitude, Julian day, climatic region, 200 m composite Tb and NDVI, and 200 m elevation, population, and land cover. We use each of these models predict the residual for all 200 m celldays (Rp_rf and Rp_gb, respectively).

2. 8 . 1 .

 81 Stage 3: increasing spatial resolution to 200 m across large urban areas In stage 3 we increase the spatial resolution of our predictions over large urban areas. We start by associating each 200 m grid cell with Tap_1km (Ta predicted in stage 2 by the final 1 km model) from the 1 km grid cell that contains the 200 m grid cell. Next, we calculate the residuals (R) for all 200 m grid cell-days with a weather station Ta observation by subtracting observed Ta from Tap_1km. The number of cell-days with a weather station observation varies by year; on average there are about 462 thousand for Tmean and 789 thousand for each of Tmin and Tmax. We use these cell-days to train a random forest and an extreme gradient boosting (XGBoost) model with the equation: R '( = 𝑓 1 T )*_,-.'( , T /'0 , NDVI '0 , Land Cover '12 , Climate ' , Elevation ' , Population ' , x ' , y ' , j I + 𝜀 '( Eq. 2 where Rij is the residual of the 1 km Ta model associated with 200 m grid cell i on day j; f designates the random forest or extreme gradient boosting function; Tap_1kmij is the 1 km Ta model prediction associated with 200 m grid cell i on day j; Tbim is the Landsat top-ofatmosphere brightness temperature of cell i for the calendar month m in which day j falls;NDVIim is the Landsat NDVI of cell i for the calendar month m in which day j falls; Land Coverily is the fraction of cell i occupied by each land cover group l in the CLC inventory year y closest to day j; Climatei is the climatic region of cell i; Elevationi is the elevation of cell i;

2. 8

 8 .2. Stage 4: improving 200 m predictions In stage 4 we improve the stage 3 predictions by ensembling. We use all 200 m grid cell-days with a weather station Ta observation to calibrate a generalized additive model (GAM) with the formula: R '( = 𝑡(x ' , y ' , ) × R *_34'( + 𝑡(x ' , y ' ) × R *_5/'( + 𝜀 '( Eq. 3 where Rij is the residual of the 1 km Ta model associated with 200 m grid cell i on day j; t(xi, yi) is a tensor product smooth of the x and y coordinates of cell i; Rp_rfij and Rp_gbij are the predicted residuals of the 1 km Ta model from the stage 3 random forest and XGBoost model, respectively, for cell i on day j; and eij is the error for cell i on day j. The GAM averages the random forest and XGBoost predicted residuals using weights that vary both by location and with the magnitude of each model's predicted residual. Finally, we add the ensemble-predicted residuals for all 200 m grid cells to Tap_1km (Ta predicted in stage 2 by the final 1 km model) to obtain daily 200 m predicted Ta (Tap_200m) across large urban areas.

Fig. 2

 2 Fig. 2 shows the spatial pattern of the daily 1 km Ta predictions of the stage 2 model on selected winter and summer days. On the cold winter day of Feb 18, 2003, predictions range from Tmin of -17 °C in parts of the Alps, the Massif Central, and the Pyrenees to Tmax of 11 °C on the Mediterranean coast. The urban heat island of Paris is faintly visible in the north center of the Tmin and Tmean maps but disappears on the Tmax map. Spatial contrasts corresponding to terrain features are well resolved, and the spatial pattern of Tmin vs. Tmean vs. Tmax varies most in the north, northeast, and southwest.

Fig. 2 .

 2 Fig. 2. Predicted 1 km Ta from the stage 2 model on selected days: Feb 18, 2003 (top row) and Aug 10, 2012 (bottom row).

Fig. 3

 3 Fig. 3 shows the spatial pattern of predicted 1 km Tmin from the stage 2 model and predicted 200 m Tmin from the stage 4 model for the Paris metropolitan area (northern France, population 12.5 million), the Toulouse metropolitan area (southwestern France, Population 1.3 million), and the Nancy metropolitan area (northeastern France, population 250,000) on the cold winter day of Feb 18, 2003. In the large city of Paris, an urban heat island is clearly visible centered over the large urban core where Tmin is about 5 °C warmer than the rural surroundings. The 200 m predictions are slightly higher than the 1 km predictions in the peripheral built-up areas and capture fine details such as the warmer Seine river and cooler parks. In the midsize city of Toulouse, the 1 km predictions capture an urban heat island over the dense city center and the suburbs to the northwest and southeast, with Tmin about 3 °C warmer than the rural surroundings. The 200 m predictions show warm Tmin in the southwestern suburbs where 1 km

Fig. 3 .

 3 Fig. 3. Predicted 1 km Tmin from the stage 2 model alone (top row) and with predicted 200 m Tmin from the stage 4 model overlaid (bottom row) on Feb 18, 2003 over the Paris, Toulouse, and Nancy metropolitan areas.

Table 1 .

 1 Satellite instruments used in this study.

	Instrument Satellite	Resolution Revisit time Overpass * Time period
	MODIS	Terra	1 km	12 hours	10:00 22:00	2000-02-02 -present
	MODIS	Aqua	1 km	12 hours	13:00 01:00	2002-07-04 -present
	TM	Landsat 5 120 m †	16 days	10:00	1984-03-01 -2011-11-18
	ETM+	Landsat 7 60 m †	16 days	10:00	1999-04-15 -present
	TIRS	Landsat 8 100 m †	16 days	10:00	2013-02-11 -present
	*Approximate local solar time; †Resampled to 30 m

Table 2 .

 2 Stage 1 model (predicting daily 1 km Ta from LST): 10-fold cross-validated performance

	across all years (2000 -2016), overall, spatial, and temporal components.
			Overall		Spatial		Temporal
		N*	R 2	RMSE MAE	R 2	RMSE MAE	R 2	RMSE MAE
	Tmin	354	0.92 1.89	1.43	0.89 1.10	0.80	0.94 1.61	1.19
	Tmean	205	0.97 1.29	0.94	0.95 0.83	0.57	0.97 1.15	0.84
	Tmax	324	0.95 1.81	1.35	0.88 1.23	0.89	0.96 1.52	1.12
	* N = mean thousands of observations used to fit each annual model

Table 3

 3 presents the 10-fold cross-validated performance of the stage 1 models across all years by calendar month and season and Table

Table 3 .

 3 Stage 1 model performance (predicting daily 1 km Ta from LST): 10-fold cross-validated performance across all years (2000 -2016), by month and season.

		Tmin			Tmean			Tmax
		R 2	RMSE MAE	R 2	RMSE MAE	R 2	RMSE MAE
	Jan	0.83 2.16	1.60	0.89 1.54	1.11	0.86 1.87	1.37
	Feb	0.84 2.03	1.51	0.91 1.37	0.99	0.89 1.74	1.28
	Mar	0.80 1.92	1.46	0.91 1.22	0.91	0.89 1.72	1.28
	Apr	0.77 1.82	1.39	0.91 1.17	0.85	0.87 1.75	1.32
	May	0.80 1.75	1.33	0.92 1.20	0.86	0.85 1.85	1.39
	Jun	0.81 1.74	1.32	0.92 1.23	0.90	0.84 1.94	1.46
	Jul	0.79 1.71	1.30	0.92 1.19	0.88	0.84 1.90	1.44
	Aug	0.78 1.77	1.35	0.92 1.18	0.88	0.87 1.89	1.43
	Sep	0.79 1.83	1.40	0.92 1.12	0.84	0.87 1.70	1.29
	Oct	0.83 1.94	1.47	0.91 1.26	0.93	0.88 1.67	1.25
	Nov	0.82 2.02	1.52	0.89 1.42	1.03	0.88 1.69	1.25
	Dec	0.82 2.17	1.61	0.86 1.69	1.21	0.84 1.94	1.39
	Winter	0.83 2.12	1.57	0.89 1.55	1.11	0.86 1.86	1.35
	Spring	0.86 1.83	1.40	0.94 1.20	0.87	0.91 1.77	1.33
	Summer	0.80 1.74	1.32	0.92 1.20	0.89	0.86 1.91	1.44
	Fall	0.87 1.92	1.46	0.95 1.26	0.92	0.93 1.69	1.27

Table 4 .

 4 Stage 1 model performance (predicting daily 1 km Ta from LST): 10-fold cross-validated performance across all years (2000 -2016), by climatic region and urban vs. rural locations.

		Tmin			Tmean			Tmax
		R 2	RMSE MAE	R 2	RMSE MAE	R 2	RMSE MAE
	Mountain	0.90 2.22	1.71	0.95 1.69	1.25	0.93 2.26	1.73
	Semi-continental	0.91 2.11	1.61	0.96 1.44	1.07	0.95 2.00	1.52
	Modified oceanic	0.94 1.53	1.16	0.98 0.98	0.73	0.98 1.33	1.01
	Transitional oceanic 0.92 1.81	1.37	0.97 1.20	0.88	0.95 1.74	1.31
	Oceanic	0.90 1.79	1.33	0.96 1.20	0.88	0.94 1.83	1.36
	Mod. Mediterranean 0.90 2.22	1.71	0.96 1.43	1.07	0.94 2.03	1.55
	Southwest basin	0.94 1.60	1.22	0.98 1.04	0.76	0.97 1.40	1.04
	Mediterranean	0.93 1.81	1.40	0.98 1.11	0.84	0.96 1.62	1.25
	Urban	0.93 1.85	1.35	0.97 1.32	0.96	0.95 1.79	1.35
	Peri-urban*	0.93 1.71	1.29	0.97 1.18	0.87	0.96 1.71	1.27
	Rural	0.92 1.90	1.44	0.97 1.30	0.94	0.95 1.82	1.36
	* non-urban locations within 5 km of a large urban area

Table 5 .

 5 Stage 4 model performance (predicting daily 200 m residuals with an ensemble): 10-fold cross-validated performance across years (2000 -2016), overall and by month and season (residual scale).

		RTmin			RTmean			RTmax
		R 2	RMSE MAE	R 2	RMSE MAE	R 2	RMSE MAE
	Overall	0.79 0.63	0.39	0.79 0.41	0.26	0.85 0.51	0.31
	Jan	0.84 0.56	0.34	0.82 0.40	0.24	0.85 0.48	0.27
	Feb	0.82 0.59	0.36	0.81 0.39	0.24	0.84 0.49	0.28
	Mar	0.80 0.63	0.39	0.79 0.40	0.26	0.83 0.50	0.30
	Apr	0.77 0.63	0.40	0.76 0.39	0.25	0.83 0.51	0.31
	May	0.77 0.60	0.37	0.76 0.38	0.24	0.84 0.51	0.31
	Jun	0.77 0.62	0.40	0.79 0.39	0.25	0.87 0.52	0.33
	Jul	0.76 0.66	0.43	0.77 0.42	0.28	0.86 0.55	0.35
	Aug	0.77 0.67	0.44	0.78 0.41	0.28	0.87 0.54	0.34
	Sep	0.77 0.69	0.46	0.75 0.42	0.29	0.84 0.54	0.34
	Oct	0.78 0.65	0.41	0.76 0.42	0.27	0.82 0.52	0.32
	Nov	0.80 0.61	0.37	0.79 0.41	0.25	0.81 0.50	0.29
	Dec	0.83 0.60	0.37	0.84 0.43	0.27	0.84 0.52	0.31
	Winter	0.83 0.58	0.36	0.83 0.41	0.25	0.84 0.50	0.28
	Spring	0.78 0.62	0.39	0.77 0.39	0.25	0.84 0.51	0.31
	Summer	0.76 0.65	0.42	0.78 0.41	0.27	0.86 0.54	0.34
	Fall	0.78 0.65	0.41	0.77 0.42	0.27	0.82 0.52	0.32

Table 6 .

 6 Stage

		RTmin			RTmean			RTmax
		R 2	RMSE MAE	R 2	RMSE MAE	R 2	RMSE MAE
	Mountain	0.83 0.67	0.42	0.83 0.46	0.30	0.88 0.58	0.36
	Semi-continental	0.81 0.66	0.42	0.79 0.43	0.28	0.86 0.55	0.34
	Modified oceanic	0.75 0.54	0.33	0.76 0.33	0.21	0.81 0.40	0.23
	Transitional oceanic	0.77 0.62	0.39	0.78 0.39	0.25	0.84 0.50	0.30
	Oceanic	0.75 0.62	0.40	0.77 0.39	0.26	0.83 0.50	0.30
	Mod. Mediterranean	0.82 0.73	0.47	0.78 0.47	0.31	0.84 0.62	0.41
	Southwest basin	0.75 0.59	0.36	0.69 0.38	0.24	0.78 0.48	0.29
	Mediterranean	0.77 0.67	0.44	0.73 0.42	0.28	0.80 0.57	0.39
	Urban	0.79 0.53	0.32	0.82 0.37	0.23	0.84 0.46	0.27
	Peri-urban*	0.76 0.58	0.36	0.78 0.37	0.24	0.83 0.47	0.28
	Rural	0.79 0.63	0.40	0.79 0.41	0.26	0.85 0.52	0.32
	* non-urban locations within 5 km of a large urban area

4 model performance (predicting daily 200 m residuals with an ensemble): 10-fold cross-validated performance across all years (2000 -2016), by climatic region and urban vs. rural locations (residual scale).

Acknowledgements

This study was funded by the Climate Health Effects In Pregnant Women And Children-A Multi-Cohort Study In France And Israel grant (CNRS PRC 2018-2020), the Fondation de France (n° 00081169), and the Effects of Urban Microclimate Variability and Global Climate Change on Heat-Related Cardiovascular Outcomes in the Semi-Arid Environment of Southern Israel grant (MOST PRC 2018-2020). Ian Hough is supported by a grant from the French National Research Agency in the framework of the "Investissements d'avenir" program (ANR-15-IDEX-02) and Ben Gurion University of the Negev. Allan C. Just is supported by NIH grants

P30ES023515 and R00ES023450. Some data processing and analyses were performed on the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA: http://www.ci-ra.org). We thank Météo France for providing data from the weather monitoring network.


Alaa Saade

Florent Krzakala

Marc Lelarge

Lenka Zdeborová

Spectral detection in the censored block model

teaching and research institutions in France or abroad, or from public or private research centers.

A. Introduction

In many inference problems, the available data can be represented on a weighted graph. Given the knowledge of the edge weights, the task is to infer latent variables carried by the nodes. Here, we shall consider the problem of recovering binary node labels from censored edge measurements [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF], [START_REF] Abbe | Conditional random fields, planted constraint satisfaction and entropy concentration[END_REF]. Specifically, given an Erdős-Rényi random graph G = (V, E) ∈ G(n, α/n) with n nodes carrying latent variables σ i = ±1, 1 ≤ i ≤ n, we draw the edge labels J ij = ±1, (ij) ∈ E from the following distribution:

P (J ij |σ i , σ j ) = (1-)1(J ij = σ i σ j )+ 1(J ij = -σ i σ j ) , (1)
where is a noise parameter. In the noiseless case = 0, we have σ i σ j = J ij and one can easily recover the communities in each connected component along a spanning tree. When = 1/2, on the other hand, the graph doesn't contain any information about the latent variables σ i , and recovery is impossible. What happens in between? The problem of exactly recovering the latent variables σ i has been studied in [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF]. It turns out that, asymptotically in the large n limit, exact recovery is shown to be possible if and only if

α > α exact = 2 log n (1 -2 ) 2 , ( 2
)
where α is the average degree of the graph. Note that the variable of an isolated vertex cannot be recovered so that the average degree has to grow at least like log n, as in the Coupon collector's problem, to ensure that the graph is connected.

We consider in this paper the case where the average degree α will remain fixed as n tends to infinity. In this setting, we cannot ask for exact recovery and we consider here a different question: is it possible to infer an assignment σi of the latent variables that is positively correlated with the planted variables σ i ? We call positively correlated an assignment σi such that the following quantity, called overlap, is strictly positive:

2 max 1 n n i=1 1(σ i = σ i ), 1 n n i=1 1(σ i = -σ i ) - 1 2 . (3)
In the limit n → ∞, this overlap vanishes for a random guess σi , and is equal to unity if the recovery is exact. We will refer to the task of finding a positively correlated assignment σi as partial recovery. This task has been shown [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF], [START_REF] Lelarge | Reconstruction in the labeled stochastic block model[END_REF] to be possible only if

α > α detect = 1 (1 -2 ) 2 . (4)
To the best of our knowledge, there is no rigorous proof that this bound is also sufficient. In [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF], the same authors also showed that belief propagation (BP) allows to saturate this bound. However, there is no rigorous analysis of BP for this problem and the fact that condition (4) is necessary and sufficient was left as a conjecture in [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF] and only the necessary part was proved in [START_REF] Lelarge | Reconstruction in the labeled stochastic block model[END_REF]. Moreover, from a practical point of view, BP requires the knowledge of the noise parameter .

In this contribution, we describe two simple spectral algorithms and we show rigorously that they are optimal, in the sense that they can perform partial recovery as soon as α > α detect . Additionally, the output of these algorithms is shown numerically to have an overlap similar to that of BP, without requiring the knowledge of the noise parameter . This closes the gap from [START_REF] Heimlicher | Community detection in the labelled stochastic block model[END_REF], [START_REF] Lelarge | Reconstruction in the labeled stochastic block model[END_REF], where spectral methods are introduced that succeed only if the connectivity is significantly larger than the threshold (4). The resulting algorithms are thus fast, trivial to implement, and asymptotically optimal.

B. Motivation and Related work

There are various interpretations and models that connect to this problem such as i) Community detection [START_REF] Abbe | Conditional random fields, planted constraint satisfaction and entropy concentration[END_REF]: we try to recover the community membership of the nodes based on noisy (or censored) observations about their relationship; ii) Correlation clustering [START_REF] Bansal | Correlation clustering[END_REF]: we try to cluster the graph G by minimizing the number of "disagreeing edges" (J ij = -1) in each cluster. These examples, and others such as synchronisation, are discussed in details in [START_REF] Abbe | Decoding binary node labels from censored edge measurements: Phase transition and efficient recovery[END_REF].

The inspiration for the present contribution comes from recent developments in the problem of detecting communities arXiv:1502.00163v2 [cs.SI] 10 Jun 2015 in the (sparse) stochastic block model. The threshold for partial recovery in the stochastic block model was conjectured in [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF] and proved in [START_REF] Mossel | Stochastic block models and reconstruction[END_REF]- [START_REF] Mossel | A proof of the block model threshold conjecture[END_REF]. Optimal spectral methods, based on the same operators as the algorithms introduced here, were proposed in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. These operators were in particular shown to be much better suited to very sparse graphs than the traditional adjacency or Laplacian operators.

Interestingly, this problem first appeared in statistical physics. Indeed, the posterior distribution corresponding to eq. (1) reads, using β 0 = 1 2 log 1-

P (σ|J) = e β0 (ij)∈E Jij σiσj Z J . (5)
This is nothing but the spin glass [START_REF] Mézard | Spin glass theory and beyond[END_REF] problem where the couplings J ij are correlated with the "planted" configuration σ i [START_REF] Abbe | Conditional random fields, planted constraint satisfaction and entropy concentration[END_REF], [START_REF] Krzakala | Hiding quiet solutions in random constraint satisfaction problems[END_REF]. Such problems can also be shown to be equivalent to spin glasses on the so-called Nishimori line [START_REF] Krzakala | Statistical physics of inference problems[END_REF], [START_REF] Nishimori | Internal energy, specific heat and correlation function of the bond-random ising model[END_REF]. With these notations, the detection condition (4) corresponds to the well-known spin glass transition [START_REF] Viana | Phase diagrams for dilute spin glasses[END_REF], [START_REF] Guerra | The high temperature region of the vianabray diluted spin glass model[END_REF] at √ α detect tanh β 0 = 1. In this spin glass context, [START_REF] Zhang | Non-backtracking operator for ising model and its application in attractor neural networks[END_REF] already conjectured that a spectral algorithm based on the nonbacktracking operator (see sec. I-A) was optimal.

C. Outline and main results

In section I, we describe two spectral algorithms that achieve the threshold (4). These algorithms are based on two linear operators: the non-backtracking operator introduced in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], and the Bethe Hessian introduced in [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. We further illustrate their properties by showing the results of numerical experiments. In section II, we list the spectral properties of the non-backtracking operator that are relevant to the present context. Finally, we discuss the properties of the Bethe Hessian and its relation with the non-backtracking operator in section III and discuss its connection with the Bethe free energy.

I. SPECTRAL ALGORITHMS

A. The non-backtracking operator

The non-backtracking operator acts on the directed edges i → j of the graph as

B i→j,k→ = J k 1(j = k)1(i = ) . (6)
It is therefore represented by a 2m × 2m matrix, where m is the number of edges in the graph. As discussed in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], [START_REF] Zhang | Non-backtracking operator for ising model and its application in attractor neural networks[END_REF] the motivation for using this operator is that it corresponds to the linear approximation of belief propagation for this problem around the so-called uninformative fixed point of BP.

Similarly to [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], one can show (see Sec. III for details) that the eigenvalues of B that are different from ±1 form the spectrum of the simpler 2n × 2n matrix

B = 0 D -1 -1 J , ( 7
)
where 1 is the n × n identity matrix, D is the diagonal matrix defined by D ii = d i , where d i is the degree of node i, and J has entries equal to the edge weights J ij . Furthermore, if

(λ = ±1, v ∈ R 2m ) is an eigenpair of B, then (λ, v ∈ R 2n ) is an eigenpair of B if v n+i = j∈∂i v j→i , ∀1 ≤ i ≤ n, (8)
λv i = (d i -1)v n+i , (9)
where ∂i and d i are the set of neighbors and the degree of node i. We will therefore favor using B . The algorithm is then as follows: given a graph with edge weights

J ij , Algorithm 1
1) build the matrix B 2) compute its leading eigenvalue λ 1 (with largest magnitude), and its corresponding eigenvector v = {v i }.

3) if λ 1 ∈ R and λ 1 > √ α
, where α is the average degree of the graph, set xi = sign(v n+i ). Otherwise, raise an error.

Theorem 1 ensures that whenever (4) holds, this algorithm outputs an assignment xi that is positively correlated with the planted latent variables x i .

B. The Bethe Hessian

Another operator closely related to the non-backtracking operator was introduced in [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. This operator, called the Bethe Hessian, is an n × n real and symmetric matrix defined as

H = (α -1)1 - √ αJ + D , (10)
where D is the diagonal matrix of vertex degrees. Based on this operator, we propose the following algorithm: given a graph with edge weights

J ij , Algorithm 2
1) build the Bethe Hessian H 2) compute its (algebraically) smallest eigenvalue λ, and its corresponding eigenvector v. 3) if λ < 0, set xi = sign(v i ). Otherwise, raise an error.

Justifications for this second algorithm, and its relation with the first one, will be provided in section III. Compared to the first algorithm, this second one is based on a smaller, symmetric matrix, which leads to improved numerical performance and stability. Additionally, in the case of more general edge weights J ij = ±1, the reduction of B to a smaller matrix B fails, and one has to work with a 2m×2m matrix. The Bethe Hessian, on the other hand, generalizes easily to arbitrary weights without any loss in scalability [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF].

C. Numerical results

Before turning to proofs, we show on figure 1 the numerical performance of our two algorithms, and compare them with the performance of belief propagation ( [3], [START_REF] Pearl | Reverend bayes on inference engines: A distributed hierarchical approach[END_REF]) which is believed to be optimal on such locally tree-like graphs in the sense that it gives, arguably, the Bayes optimal value of the overlap asymptotically. As shown in section II, both algorithms 1 and 2 are able to achieve partial recovery as soon as α > α detect , and their overlap is similar to that of BP, though of course strictly smaller. Note again that BP requires the knowledge of while the two spectral algorithms described here do not, are The overlap for BP is estimated asymptotically using the standard method of population dynamics (see for instance [START_REF] Mezard | Information, physics, and computation[END_REF]), with a population of size 10 4 . All three methods output a positively correlated assignment as soon as α > α detect . Spectral algorithms 1 and 2 have an overlap similar to that of BP, with the same phase transition, while being simpler and not requiring the knowledge of the parameter .

trivial to implement, run faster, and avoid the potential nonconvergence problem of belief propagation while remaining asymptotically optimal in detecting the hidden assignment. We also observe, empirically, that the overlap given by the Bethe Hessian seems to be always superior to the one provided by the non-backtracking operator.

II. SPECTRAL PROPERTIES OF THE NON-BACKTRACKING

OPERATOR

In this section, we state results concerning the spectrum of B and show that algorithm 1 outputs an assignment σi that is positively correlated with the planted one, whenever (4) holds.

As already noticed in previous work for the case of an unweighted random graph [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], [START_REF] Saade | Spectral density of the non-backtracking operator on random graphs[END_REF], the superior performance of the non-backtracking operator B is due to the particular shape of its spectrum. In the case of the stochastic block model [START_REF] Holland | Stochastic blockmodels: First steps[END_REF], it decomposes into a bulk of uninformative eigenvalues contained in a disk of radius √ α in the complex plane, and a few real and informative eigenvalues outside of the disk. This observation was recently proven in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF], in the case of 2 communities.

The following theorem generalizes this previous result to the present setting and is the main result of this paper.

Theorem 1: Given an Erdős-Rényi random graph with average degree α, variables assigned to vertices σ i = ±1 uniformly at random independently from the graph and where the edges carry weights sampled from (1), we denote by B the non-backtracking operator defined by [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF]. and by

|λ 1 | ≥ |λ 2 | ≥ • • • ≥ |λ 2m |
the eigenvalues of B in order of decreasing magnitude. Then, with probability tending to 1 as n → ∞, we have:

(i) if α < α detect then |λ 1 | ≤ √ α + o(1). (ii) if α > α detect , then λ 1 ∈ R, λ 1 = α(1-2 )+o(1) > √ α, and |λ 2 | ≤ √ α + o(1).
Additionally, denoting v the eigenvector associated with λ 1 , the following assignment is positively correlated with the planted variables σ i :

σi = sign   j∈∂i v j→i   .
This theorem is illustrated on Fig. 2.

It is then straightforward to show the following:

Corollary 1: The assignment output by Algo. 1 is positively correlated with the planted variables σ i if and only if

α > α detect . ( 11
)
We now give a brief sketch of proof for our Theorem 1. The proof relies heavily on the techniques developed in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]. We try to use notation consistent with [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]: E is the set of oriented edges and for any e = u → v = (u, v) ∈ E, we set e 1 = u, e 2 = v and e -1 = (v, u). For a matrice M , its transpose is denoted by M * . We start with a simple observation: if t is the vector in R E defined by t e = σ e2 and is the Hadamard product, i.e. (t x) e = σ e2 x e , then we have

Bx = λx ⇔ B(t x) = λ(t x), (12)
with B defined by Bef = B ef σ f1 σ f2 . In particular, B an B have the same spectrum and there is a trivial relation between their eigenvectors. It will be easier to work with B so to lighten the notation, we will denote (in this section):

B ef = 1(e 2 = f 1 )1(e 1 = f 2 )P f ,
where P f = σ f1 J f σ f2 . Note that the random variables P f are now i.i.d. with P(P f = 1) = 1 -P(P f = -1) = 1 -. With this formulation, the problem is said in statistical physics to be "on the Nishimori line" [START_REF] Krzakala | Statistical physics of inference problems[END_REF], [START_REF] Nishimori | Internal energy, specific heat and correlation function of the bond-random ising model[END_REF].

For the case (1 -2 ) 2 α < 1, the proof is relatively easy. Indeed, from [START_REF] Lelarge | Reconstruction in the labeled stochastic block model[END_REF], we know that our setting is contiguous to the setting with = 1/2. In this case, the random variable P i,j are centered and a version of the trace method will allow to upper bound the spectral radius of B. Note however, that one needs to condition on the graph to be -tangle-free, i.e. such that every neighborhood of radius contains at most one cycle in order to apply the first moment method.

We now consider the case (1 -2 ) 2 α > 1 and denote by P the linear mapping on R E defined by (P x) e = P e x e -1 (i.e. the matrix associated to P is P ef = P e 1(f = e -1 )). Note that P * = P and since P 2 e = 1, P is an involution so that P is an orthogonal matrix. A simple computation shows that B k P = P B * k , hence B k P is a symmetric matrix. This symmetry corresponds to the oriented path symmetry in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF] and will be crucial to our analysis.

We also define α = (1 -2 )α and χ ∈ R E with χ e = 1 for all e ∈ E. The proof strategy is then similar to Section 5 in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]. Consider a sequence ∼ κ log α n for some small positive κ. Let

ϕ = B χ B χ , θ = B P ϕ , ζ = B P ϕ θ .
If R = B -θζP ϕ * and we can prove that R is small in comparison with θ, then we can use a theorem on perturbation of eigenvalues and eigenvectors adapted from the Bauer-Fike theorem (see Section 4 in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]) saying that B should have an eigenvalue close to θ.

More precisely, for y ∈ R E with y = 1, write y = sP ϕ + x with x ∈ (P ϕ) ⊥ and s ∈ R. Then, we find

Ry = B x + s(B P ϕ -θζ) ≤ sup x: x,P ϕ =0, x =1 B x .
This last quantity can be shown to be upper bounded by (log n) c α /2 similarly as in Proposition 12 in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]. Moreover, we can also show that w.h.p.

ζ, P ϕ ≥ c 0 , c 0 α ≤ θ ≤ c 1 α . (13)
These bounds allow to show that B has an eigenvalue

|λ 1 - α| = O(1/ ) and that |λ 2 | ≤ √ α + o(1). Note that θ = B B * P χ B χ
, so that we need to compute quantities of the type B χ . We now explain the main ideas to compute these quantities. First note that, (B χ) e depends only on the ball of radius around the edge e. For not too large, this neighborhood can be coupled with a Galton-Watson branching process with offspring distribution Poi(α). It is then natural to consider this Poisson Galton-Watson branching process with i.i.d. weights P u,v ∈ {±1} on its edges with mean 1 -2 . For u in the tree, we denote by |u| its generation and by Y (u) = t s=1 P γs,γs+1 where γ = (γ 1 , . . . , γ ) is the unique path between the root o = γ 1 and u = γ t . Then (B χ) e is well approximated by:

Z = |u|= Y (u).
It is easy to see that X t = Zt αt is a martingale (with respect to the natural filtration) with zero mean. Moreover we have

E Z 2 t = E   u,v:|u|=|v|=t Y (u)Y (v)   = t i=0 α t-i (1 -2 ) 2i α 2i = O α2t ,
where the last equality is valid only if (1 -2 ) 2 α > 1. So in this case, we have E X 2 t = O(1) and the martingale X t converges a.s. and in L 2 to a limiting random variable X(∞) with mean one. Following the argument as in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF], this reasoning leads to [START_REF] Krzakala | Hiding quiet solutions in random constraint satisfaction problems[END_REF].

We now consider the eigenvector associated with λ 1 . It follows from Bauer-Fike theorem (see Section 4 in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]) that the eigenvector x associated to λ 1 is asymptotically aligned with B B * P χ B B * P χ . Thanks to the coupling with the branching process, we can prove that B B * P χ ≈ α2 and moreover, we have for e ∈ E,

(B B * P χ) e α2 ≈ α α(1 -2 ) 2 -1 X(∞), (14)
where X(∞) is the limit of the martingale defined above and has mean one. We can now translate this result to the eigenvector of the original non-backtracking operator thanks to [START_REF] Mézard | Spin glass theory and beyond[END_REF]: v e = σ e2 x e where x e is approximated by [START_REF] Krzakala | Statistical physics of inference problems[END_REF]. In particular, we see that e,e2=v v e is correlated with σ v .

III. FROM THE NON-BACKTRACKING OPERATOR TO THE BETHE HESSIAN

In this section, we relate the spectra of H, B and B by generalizing some properties discussed in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. (λ = ±1, v ∈ R 2m ) being an eigenpair of B, we define

v i = j∈∂i v j→i , ∀1 ≤ i ≤ n . (15)
Since λv i→j = k∈∂i\j J ki v k→i it follows that λv i→j = v i -J ij v j→i . Closing the equation on the single site elements v i thus leads to

v i 1 + k∈∂i J 2 ij λ -J 2 ij -λ k∈∂i J ij λ -J 2 ij v k = 0 . (16)
For convenience, we now define the matrix:

H(X) = (X 2 -1)1 -XJ + D (17)
Note in particular that the Bethe Hessian reads H = H( √ α). Given that the values of J ij are ±1, all eigenvalues of B different from ±1 thus must satisfies the following generalization of the Ihara-Bass formula [START_REF] Bass | The ihara-selberg zeta function of a tree lattice[END_REF] :

det (λ 2 -1)1 -λJ + D = det H(λ) = 0 . (18)
To solve [START_REF] Viana | Phase diagrams for dilute spin glasses[END_REF] one needs to find an eigenvector v of H(λ) with a zero eigenvalue. This is a quadratic eigenproblem, which can be turned into a linear one by introducing the matrix B of Algo. 1. Indeed, if λ ∈ R is an eigenvalue of B with eigenvector v , then it follows that v := {v i } n+1≤i≤2n is an eigenvector of H(λ) with eigenvalue 0, so that λ is an eigenvalue of B as well (at least if λ = ±1), justifying eq. [START_REF] Massoulie | Community detection thresholds and the weak ramanujan property[END_REF][START_REF] Mossel | A proof of the block model threshold conjecture[END_REF]. Note that since we are interested in values of λ > 1 (since λ > α and we need α > 1 from (4)), the limitation of looking at λ = ±1 is irrelevant.

Finally, following [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], we can relate the spectra of B and H by the following argument. For X large enough, H(X) is positive definite. Then as X decreases, H(X) will gain a new negative eigenvalue whenever X becomes equal to an eigenvalue of B. This justifies the following corollary:

Corollary 2: if the conditions of Theorem 1 apply, then H = H( √ α) has a unique negative eigenvalue if α > α detect , and none otherwise.

Strictly speaking, if we denote by λ 1 the leading eigenvalue of B, we have only shown that the eigenvector with eigenvalue 0 of H(λ 1 ) is positively correlated with the planted variables if α > α detect . However, we observe numerically (see figure 1) that the eigenvector with negative eigenvalue of H is also positively correlated, and in fact gives a slightly better overlap. This point will have to be clarified in future work.

It is worth noting the Bethe Hessian is also related to the belief propagation algorithm. [START_REF] Yedidia | Bethe free energy, kikuchi approximations, and belief propagation algorithms[END_REF] showed that the fixed points of the BP recursion are stationary points of the socalled Bethe free energy. Direct optimization of the Bethe free energy has then been proposed as an alternative to BP [START_REF] Welling | Belief optimization for binary networks: A stable alternative to loopy belief propagation[END_REF]. In this context, [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF] showed that the so-called paramagnetic fixed point (corresponding to an uninformative assignment) is a local minimum of the Bethe free energy if and only if H is positive definite. Algo. 2 can therefore be seen as a spectral relaxation of the direct optimization of the Bethe free energy. In the end, both approaches are indeed deeply related to BP.

IV. CONCLUSION

We have considered the problem of partially recovering binary variables from the observation of censored edge weights, and described two optimal spectral algorithms for this task that can provably perform partial recovery as soon as it is information theoretically possible to do so. Remarkably, these algorithms do not require the knowledge of the noise parameter and perform almost as well as belief propagation, which is expected (but not proved) to be Bayes optimal for this problem. This allows to close the gap from previous works, both algorithmically, by providing optimal spectral algorithms, and theoretically, by proving that the transition (4) is a necessary and sufficient condition for partial recovery.

Fig. 1 .

 1 Fig.1. Overlap as a function of α: comparison between algorithm 1 (based on the non-backtracking operator B), algorithm 2 (based on the Bethe Hessian H), and belief propagation (BP). The noise parameter is fixed to 0.25 (corresponding to α detect = 4), and we vary α. The overlap for B and H is averaged over 20 graphs of size n = 10 5 . The overlap for BP is estimated asymptotically using the standard method of population dynamics (see for instance[START_REF] Mezard | Information, physics, and computation[END_REF]), with a population of size 10 4 . All three methods output a positively correlated assignment as soon as α > α detect . Spectral algorithms 1 and 2 have an overlap similar to that of BP, with the same phase transition, while being simpler and not requiring the knowledge of the parameter .

Fig. 2 .

 2 Fig.2. Spectrum of the non-backtracking matrix in the complex plane for a problem generated with = 0.25, n = 2000. We used α = 3 (left side) and α = 8 (right side), to be compared with α detect = 4. Each point represents an eigenvalue. In both cases, the bulk of the spectrum is confined in a circle of radius √ α. However, when α > α detect , a single isolated eigenvalue appears out of the bulk at (1 -2 )α (see the arrow on the right plot) and the corresponding eigenvector is correlated with the planted assignement.

ACKNOWLEDGMENT

This work has been supported by the ERC under the European Union's FP7 Grant Agreement 307087-SPARCS and by the French Agence Nationale de la Recherche under reference ANR-11-JS02-005-01 (GAP project).


Shotaro Uehara

Makiko Shimizu

Karen Plé

Sylvain Routier

Nao Yoneda

Yuichiro Higuchi

Hiroshi Suemizu

Hiroshi Yamazaki
email: hyamazak@ac.shoyaku.ac.jp

Probe drug T-1032 N-oxygenation mediated by cytochrome P450 3A5 in human hepatocytes in vitro and in humanized-liver mice in vivo

Keywords: CYP3A5.1, CYP3A5*1, hepatocytes, N-oxide, T-5

Introduction

Approximately half of new medicines recently approved by the US Food and Drug Administration are reportedly metabolized by cytochrome P450 3A4 (CYP3A4) and 3A5 [START_REF] Bhutani | FDA Approved Drugs from 2015-June 2020: A Perspective[END_REF].

CYP3A4 and 3A5 have significant substrate overlap [START_REF] Niwa | Drug interactions between nine antifungal agents and drugs metabolized by human cytochromes P450[END_REF][START_REF] Yamazaki | Potential impact of cytochrome P450 3A5 in human liver on drug interactions with triazoles[END_REF], but there is currently no way to selectively monitor the activities of these two CYP3A enzymes; instead, the cumulative activity is often attributing to CYP3A4. Whereas CYP3A4 expression is ubiquitous, CYP3A5 is expressed variably in different populations, resulting in large individual differences in CYP3A5 expression levels [START_REF] Yamaori | CYP3A5 Contributes significantly to CYP3A-mediated drug oxidations in liver microsomes from Japanese subjects[END_REF].

T-1032 [methyl (2-(4-aminophenyl)-1,2-dihydro-1-oxo-7-(2-pyridinylmethoxy)-4-(3,4,5trimethoxyphenyl)-3-isoquinoline carboxylate sulfate] was synthesized as a selective phosphodiesterase 5 inhibitor [START_REF] Ukita | Novel, potent, and selective phosphodiesterase 5 inhibitors: synthesis and biological activities of a series of 4-aryl-1-isoquinolinone derivatives[END_REF]. It has been suggested that the oxidation of T-1032 (also known as T-5) by polymorphic CYP3A5 could yield an N-oxide metabolite that may be more than 100-fold selective compared with P450 3A4 in human liver microsomal P450 forms in in vitro systems [6]. Against this background, we set out to investigate whether T-1032 is a suitable CYP3A5 probe in human hepatocytes in vitro and in humanized-liver mice in vivo.

The identification of suitable CYP3A5 probes could contribute to drug metabolism and drug interaction studies. We report herein the significant but limited contributions of polymorphic CYP3A5 to the oxidation of T-1032 in in vivo systems. transplanted into NOG-TKm30 mice [START_REF] Uehara | An improved TK-NOG mouse as a novel platform for humanized liver that overcomes limitations in both male and female animals[END_REF] to create humanized-liver mice. Cryopreserved hepatocytes from 12 humanized-liver mice transplanted with human hepatocytes from the above four donors (Hu-Liver cells) were prepared using the collagenase perfusion method described previously [START_REF] Higuchi | Functional polymer-dependent 3D culture accelerates the differentiation of HepaRG cells into mature hepatocytes[END_REF][START_REF] Uehara | Expression and inducibility of cytochrome P450s in human hepatocytes isolated from chimeric mice with humanised livers[END_REF]. The percentages of human leukocyte antigen-positive cells in hepatocytes from these humanized liver mice were all >95%, as determined by fluorescenceactivated cell sorting using a monoclonal anti-human leukocyte antigen antibody (clone G46-2.6; BD Biosciences, Franklin Lakes, NJ, USA) on a BD FACSAria cell sorter (BD Biosciences). This study was reviewed and approved by the Institutional Animal Care and Use

Materials and methods

In vitro analyses

Committees of the Central Institute for Experimental Animals (Permit Number: 20060 A).

Cell suspension metabolic assays using cryopreserved hepatocytes from one human donor and fresh hepatocytes (1 × 10 5 cells/tube) from the 12 humanized-liver mice were incubated for 60 min at 37°C in an incubator on an orbital shaker. The reaction mixture (0.2 mL) contained T-1032 (10 µM) and Williams' Medium E without phenol red (Gibco, Grand Island, NY, USA) [START_REF] Uehara | Methyl-hydroxylation and subsequent oxidation to produce carboxylic acid is the major metabolic pathway of tolbutamide in chimeric TK-NOG mice transplanted with human hepatocytes[END_REF]. Reactions were started by adding the substrate and stopped by adding an equal volume of ice-cold acetonitrile containing the deuterium-labelled internal standards.

Precipitates were then removed by centrifugation (10 min at 20,000 × g, 4°C). Supernatants (2.0 μL) were analysed using liquid chromatography-tandem mass spectrometry.

In vivo animal experiments

Humanized-liver NOG-TKm30 mice (4-5 males per group, weight ~20-30 g, Central Institute for Experimental Animals, Kawasaki, Japan [START_REF] Uehara | An improved TK-NOG mouse as a novel platform for humanized liver that overcomes limitations in both male and female animals[END_REF]) were used in this in vivo study. The human hepatocyte replacement indices of the host livers were all >80%, as estimated using a human albumin ELISA quantitation kit (Bethyl Laboratories, Montgomery, TX, USA). Plasma samples were collected 0.25-24 h after oral doses of T-1032 (1.0 mg/kg). The pharmacokinetic profiles of T-1032 and its N-oxide were determined between 0.25 and 24 h using liquid chromatography and tandem mass spectrometry.

A Nexera ultra-high-performance liquid chromatography system (Shimadzu, Kyoto, Japan) Statistical analyses of the pharmacokinetic parameters of T-1032 and its N-oxide in humanized-liver mice were performed using two-way analysis of variance and unpaired t-tests (Prism, GraphPad Prism Software, La Jolla, CA, USA).

Results and discussion

To evaluate in vitro the roles of polymorphic CYP3A5 in the hepatic N-oxidation of T-1032 in humans and humanized-liver mice, T-1032 N-oxide formation rates were investigated after incubation with human hepatocytes. T-1032 N-oxide was produced by all the hepatocytes tested (Fig. 1). The hepatocytes from one human donor genotyped as CYP3A5*3/*3 (poor expressers) showed significantly lower T-1032 N-oxidation rates than those from donors harboring the CYP3A5*1 allele (extensive expressers, p < 0.01).

The in vivo pharmacokinetics of T-1032 and its N-oxygenated metabolite were investigated after administrations in humanized-liver mice. Orally administered T-1032 (1.0 mg/kg) was oxidized to its N-oxide in humanized-liver mice and was detected in plasma (Fig. 2A and2B).

T- 1032 (

 1032 CAS no. 212500-03-3), its N-oxide, T-1032-d3, and T-1032 N-oxide-d3 were obtained from Biopredic International, Saint Grégoire, France. Cryopreserved hepatocytes from four human donors (a 2-year-old Caucasian boy [donor A], a 35-year-old Asian man [donor B], a 30-year-old African American woman [donor C], and a 12-year-old Caucasian girl [donor D], genotyped as CYP3A5*1/*3, 3A5*1/*3, 3A5*1/*7, and 3A5*3/*3, respectively) were obtained from Biopredic International or Lonza Walkersville, Walkersville, MD, USA, and were

  equipped with an API QTRAP 5500 triple quadrupole mass spectrometer (AB Sciex, Framingham, MA, USA) and a YMC-Triart C18 column (100 mm × 3.0 mm i.d., particle size 3.0 μm; YMC, Kyoto, Japan) was used. For internal standards, T-1032-d3 and T-1032 N-oxide-d3 were used. Levels of T-1032 and its N-oxide in mouse plasma along with levels of T-1032-d3 and T-1032 N-oxide-d3 were quantified using a gradient elution program followed by a tandem mass spectrometry system with multiple reaction monitoring and positive ion mode electrospray ionization using the following transitions: m/z 568→476, m/z 584→476, m/z 571→479, and m/z 587→479, respectively. The mobile phase had a flow rate of 0.3 mL/min with a multistep gradient elution program using 0.1% formic acid in water (A) and methanol (B) as follows: 0-12 min, 20%-80% B; 12.1-15 min, 100% B; and 15.1-18 min, 20% B. The column temperature was maintained at 40°C throughout the run. The retention times for authentic T-1032 and its N-oxide were 10.3 and 9.1 min, respectively, under the present conditions. The pharmacokinetic parameters for T-1032 and its N-oxide were derived from the respective plasma concentration-time curves (0.25-24 h after oral doses) using trapezoidal linear interpolation methods with uniform weighting and non-compartmental analysis in Phoenix 8.2 (Pharsight, Mountain View, CA, USA).

  The area under the concentration-time curve of T-1032 N-oxide was one-tenth that of T-1032 in humanized-liver mice. Concentrations of T-1032 N-oxide in plasma from humanized-liver mice with CYP3A5*1/*7 human hepatocytes were significantly higher than those in plasma from humanized-liver mice with CYP3A5*3/*3 hepatocytes by two-way analysis of variance (p < 0.01, Fig.2B). The nonparametric pharmacokinetic parameters for T-1032 and its N-oxygenated metabolite in individual humanized-liver mice in vivo are summarized in Fig. 2C and 2D. The mean maximum concentrations (Cmax) and standard deviations in plasma from humanized-liver mice with CYP3A5*3/*3 and CYP3A5*1/*7 hepatocytes were calculated to be 0.069 ± 0.024 nmol T-1032/mL and 0.005 ± 0.001 nmol T-1032 N-oxide/mL and 0.092 ± 0.053 nmol T-1032/mL and 0.011 ± 0.005 nmol T-1032 N-oxide/mL, respectively. The mean areas under the concentration-time curves (AUC) in plasma form humanized-liver mice with CYP3A5*3/*3 and CYP3A5*1/*7 hepatocytes, respectively, were calculated to be 0.41 ± 0.13 nmol T-1032 h/mL and 0.042 ± 0.018 nmol T-1032 N-oxide h/mL and 0.45 ± 0.11 nmol T-1032 h/mL and 0.066 ± 0.023 nmol T-1032 N-oxide h/mL. After oral administration of T-1032, the Cmax values of T-1032 N-oxide in plasma from humanized-liver mice with CYP3A5*1/*7 hepatocytes were twice (p < 0.05) those from humanized-liver mice with CYP3A5*3/*3 hepatocytes (Fig. 2D). In this study, the pharmacokinetic profiles of T-1032 and its N-oxide were determined up to 24 h post-administration in humanized-liver mice. The ratio of the AUC of T-1032 N-oxide to the sum of the AUCs of T-1032 and its N-oxide were 13% (0.066 nmol T-1032 N-oxide h/mL compared with a total of 0.519 ± 0.133 nmol h/mL) in plasma from humanized-liver mice with CYP3A5*1/*7 hepatocytes. In contrast, the ratio of the AUC of T-1032 N-oxide to the sum of the AUCs of T-1032 and its N-oxide were 9.2% (0.042 nmol T-1032 N-oxide h/mL compared with a total of 0.452 ± 0.144 nmol h/mL) in plasma from humanized-liver mice with CYP3A5*3/*3 hepatocytes (poor expressers). The data shown in Figs. 2B and 2D imply significant but little effects of wild-type CYP3A5.1 on the total CYP3A activities in in vivo systems. In vitro T-1032 deletion mediated mainly by recombinant CYP3A5 and minorly by CYP2C8, 2C9, and 3A4 has been reported, probably forming the O-dealkylated metabolite

Figure legends Figure 1 .

 legends1 Figure legends

Figure 2 .

 2 Figure 2. Plasma concentrations (A, B) and pharmacokinetic parameters (C, D) for T-

Acknowledgments

We are grateful to Drs. Christophe Chesne at Biopredic and Pierre Daligaux at Starlight for their kind support of this study. The authors thank Hiroaki Kato, Yasuhiko Ando, Takaya Honma, and Norie Murayama for their technical help and David Smallbones for copyediting a draft of this article. This work was supported partly by the Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research 20K06463 and 20K07164. K.P and S.R. gratefully acknowledge funding from the Centre Val de Loire Region, France (APR-IR FLUPO), as well as the RTR Motivhealth (2019-00131403), the Labex SYNORG (ANR-11-LABX-0029) and IRON (ANR-11-LABX-0018-01), and the Feder Techsab (FEDER-FSE

Conflicts of interest

All authors declare that they have no conflicts of interest.

2014-2020-EX011313).

Author contributions H.S. and H.Y. designed the study; K.P. and S.R. optimized chemical synthesis of T-1032 and its N-oxide as well as the design and synthesis of all deuterated analogues; S.U., M.S., N.Y., Y.H., and H.S. performed the experiments; and M.S. and H.Y. modeled the pharmacokinetics.

All co-authors performed data analysis. S.U. and H.Y. wrote the manuscript.