Kamel Lahssini
email: kamel.lahssini@inrae.fr

Florian Teste
email: fteste96@protonmail.com

Karun Reuel Dayal
email: karun.dayal@inrae.fr

Sylvie Durrieu
email: sylvie.durrieu@inrae.fr

Member, IEEE Dino Ienco
email: dino.ienco@inrae.fr

Jean-Matthieu Monnet
email: jean-matthieu.monnet@inrae.fr

Sylvie Durrie

Combining LiDAR Metrics and Sentinel-2 Imagery to Estimate Basal Area and Wood Volume in Complex Forest Environment via Neural Networks

Keywords: Convolutional neural networks (CNNs), forest monitoring, multiscale remote sensing, multisensor data fusion, structure and biophysical variables estimation. I

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

terrestrial carbon stock is contained in these natural carbon sinks [START_REF] Pan | A large and persistent carbon sink in the world's forests[END_REF]. Moreover, they are home to a rich, but also very fragile, animal and plant biodiversity [START_REF] Simberloff | The role of science in the preservation of forest biodiversity[END_REF]. Their sustainable management and conservation is therefore a crucial challenge in the current context of global warming.

To address these issues and provide stakeholders with valuable inputs to support the decision-making process, the effective exploitation of available and up-to-date earth observation data is a major challenge. Nowadays, ongoing space missions and airborne acquisition campaigns allow to collect an increasing amount of remote-sensing data [START_REF] Boyd | Satellite remote sensing of forest resources: Three decades of research development[END_REF] that provide complementary information acquired via various modalities and at different spatial and temporal resolutions [START_REF] Janssen | Principles of Remote Sensing: An Introductory Textbook[END_REF].

These remote-sensing data supply valuable information to monitor and characterize forest ecosystems [START_REF] Defries | Land use change around protected areas: Management to balance human needs and ecological function[END_REF]- [START_REF] Goff | Sustainable management of old-growth boreal forests: Myths, possible solutions and challenges[END_REF].

For instance, light detection and ranging (LiDAR) technology is well suited for the characterization of structure and biophysical variables like the amount of wood resource and biomass [START_REF] Franklin | Ecological characteristics of old-growth Douglas-fir forests[END_REF], [START_REF] Chave | Tree allometry and improved estimation of carbon stocks and balance in tropical forests[END_REF], due to the fact that the laser signal can penetrate through forest canopies. Airborne laser scanning (ALS) systems are LiDAR-based systems which provide 3-D multiecho point clouds describing the structure of the vegetation from the top of the canopy to the ground. The 3-D spatial distribution of the points resulting from the laser-environment interaction, thus, can be directly linked with the 3-D spatial distribution of the vegetation [START_REF] Magnussen | Derivations of stand heights from airborne laser scanner data with canopy-based quantile estimators[END_REF]. For regional-scale monitoring, the analysis of ALS data mainly relies on area-based approaches (ABA) involving statistical analysis of the spatial distribution of the points at stand level, i.e., around 300 to 700 m 2 [START_REF] Hyyppä | Review of methods of small-footprint airborne laser scanning for extracting forest inventory data in boreal forests[END_REF]. The joint use of ground truth (GT) information at plot level and metrics derived from the ALS point clouds enables to develop empirical models to estimate structure and biophysical forest variables. Models are then applied to the whole area of study to produce forest resource maps [START_REF] Bouvier | Generalizing predictive models of forest inventory attributes using an area-based approach with airborne LiDAR data[END_REF], [START_REF] Véga | On the interest of penetration depth, canopy area and volume metrics to improve Lidar-based models of forest parameters[END_REF]. ABA models have exhibited robustness and effectiveness in predicting inventory variables for a wide range of homogeneous forests [START_REF] Nasset | Effects of different flying altitudes on biophysical stand properties estimated from canopy height and density measured with a small-footprint airborne scanning laser[END_REF]. However, they exhibit limitations when implemented in more complex and heterogeneous environments requiring specific calibration to identify the most relevant variables that best describe a given study area through ALS data [START_REF] Vincent | Accuracy of small footprint airborne LiDAR in its predictions of tropical moist forest stand structure[END_REF].

When analyzing ALS point clouds, research studies [START_REF] Zolkos | A meta-analysis of terrestrial aboveground biomass estimation using lidar remote sensing[END_REF] indicate that topography, amongst other factors, greatly influences ALS-based models' performances. This is mainly due to the fact that ALS point clouds are usually normalized in height to take into account ground elevation. The normalization process induces a distortion of the tree structure and directly affects the data distribution within each tree, influencing the metrics derived from the point clouds [START_REF] Vega | PTrees: A point-based approach to forest tree extraction from lidar data[END_REF].

Sentinel-2 optical images provide radiometric information that can be useful for the characterization of stand composition [START_REF] Karasiak | Mapping tree species of forests in southwest France using Sentinel-2 image time series[END_REF], [START_REF] Grabska | Forest stand species mapping using the Sentinel-2 time series[END_REF] and, thus, bring complementary information to ALS geometric data. The authors of [START_REF] Wittke | Comparison of two-dimensional multitemporal Sentinel-2 data with three-dimensional remote sensing data sources for forest inventory parameter estimation over a boreal forest[END_REF] investigated the usability of Sentinel-2 2-D images in forest inventories against other remote-sensing source-providing 2.5-or 3-D-data-like ALS, elevation model induced by high-resolution optical satellite WorldView-2 images, and syntheticaperture radar stereo data from TerraSAR-X. The results confirmed that the higher spatial resolution input data, in the case of 3-D products, correlate with more accurate forest inventory parameter predictions, which is in line with other results presented in literature. However, the authors also highlight that, despite not having outstanding performances, Sentinel-2 imagery has the advantage, compared to other remote-sensing data sources, to being free of charge and regularly updated. The authors of [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF] evaluated a particular machine learning algorithm Bayesian additive regression trees (BART) to analyze Sentinel-2 images and topographic variables to estimate the forest stand characteristics. The results indicate that the combination of radiometric and topographic variables (derived from PALSAR data) improved the estimation of the forest attributes compared to the use of only Sentinel-2 information.

Several studies demonstrated the value of combining LiDAR data with radiometric information from high spatial resolution multispectral satellite sensors regarding the characterization of both complex forest environments [START_REF] Monnet | Estimation de paramétres forestiers par Lidar aéroporté et imagerie satellitaire RapidEye: étude de sensibilité[END_REF], [START_REF] Wu | Estimating forest and woodland aboveground biomass using active and passive remote sensing[END_REF] and urban vegetation [START_REF] Zhang | Assessing of urban vegetation biomass in combination with LiDAR and high-resolution remote sensing images[END_REF]. However, enriching ALS-derived metrics with additional radiometric information from optical images and topography descriptors is challenging in commonly used ABA models based on linear regression. Indeed, such models do not allow to fully exploit the complex interplay between all the available input modalities since all the information is treated equally.

In recent years, deep learning (DL) approaches [START_REF] Lecun | Deep learning[END_REF] are getting more and more attention in the field of remote sensing [START_REF] Reichstein | Deep learning and process understanding for datadriven Earth system science[END_REF], since they are demonstrating compelling performances in the analysis of multimodal (i.e., multisource, multitemporal,and multiscale) earth-observation data [START_REF] Hong | More diverse means better: Multimodal deep learning meets remote-sensing imagery classification[END_REF]. Neural networks, through the stacking of nonlinear functions, allow representation and modeling of complex relationships between a set of input and output variables in an end-to-end manner. In a multisource remote-sensing scenario, two main strategies, namely 1) early fusion and 2) late fusion, are mainly adopted in order to combine heterogeneous and complementary information [START_REF] Hong | More diverse means better: Multimodal deep learning meets remote-sensing imagery classification[END_REF]. While the former strategy firstly combines data together and then feeds a standard neural network architecture, the latter one firstly analyzes each source by means of a dedicated branch (encoder) and then, terminates the processing of the data, combining the intermediate per-branch information via additional neural network layers. In both strategies, the set of network parameters are optimized end-to-end. The majority of the proposed approaches were introduced to cope with the issue of multimodal remote-sensing analysis for classification purposes [START_REF] Interdonato | DuPLO: A dual view point deep learning architecture for time series classification[END_REF]- [START_REF] Benedetti | M 3 Fusion: A. deep learning architecture for multiscale multimodal multitemporal satellite data fusion[END_REF]. Regarding the specific combination of spectral optical and Li-DAR remote-sensing data, several works exist that deal with the general task of land cover mapping by employing both earlyand late-fusion strategy to combine hyperspectral and LiDAR data [START_REF] Ghamisi | Hyperspectral and LiDAR data fusion using extinction profiles and deep convolutional neural network[END_REF]- [START_REF] Hang | Classification of hyperspectral and LiDAR data using coupled CNNs[END_REF], multispectral (three or four spectral bands) with LiDAR information [START_REF] Sun | Extracting building boundaries from high resolution optical images and lidar data by integrating the convolutional neural network and the active contour model[END_REF], or hyperspectral, multispectral, and LiDAR data together [START_REF] Bigdeli | An ensemble deep learning method as data fusion system for remote sensing multisensor classification[END_REF].

However, there is a lack of methodological investigation to leverage the potential of DL-based strategy for regression tasks in the context of multimodal remote sensing analysis. Only recently, the computer science and signal processing communities have started to rigorously investigate the potential of such strategies for monosource regression tasks [START_REF] Lathuilière | A comprehensive analysis of deep regression[END_REF], while some attempts related to the characterization of forest properties are emerging in the remote-sensing community [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF]. In [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF], the authors proposed a stacked sparse autoencoder network from LiDAR metrics and optical indices (from Landsat 8 imagery) with the aim to estimate forest aboveground biomass. In this approach, all the information (LiDAR metrics, optical indices, and combined optical-LiDAR indices) are fed to a fully-connected neural network model, following an early-fusion strategy, in order to estimate the biomass quantity.

To cope with the estimation of structure and biophysical forest variables from multimodal remote-sensing data, we propose a DL-based model capable to fully exploit the interplay among multimodal information coming from Sentinel-2 images and derived from ALS point clouds. Conversely to previous attempts in the literature, the proposed framework, named multimodal forest variables estimation based on DL framework (MMFVE), implements łate-fusion strategy deploying a two-branch DL architecture to deal with multimodal remote-sensing data: Information derived from ALS point clouds on one hand and time series of Sentinel-2 optical imagery on the other hand. We focus on the estimation of two forest variables, 1) total volume and 2) basal area, in the study area of Massif des Bauges Natural Regional Park (France). The total volume of a tree is defined as the volume of all the wood contained in it (stem and branches). When coupled with additional information, it is useful for the computation of biomass. The basal area of a tree is the cross-sectional area of the tree at breast height (usually defined as 1.3 m above ground). When computed for all the trees in a plot (m 2 /ha), it is the measure of the cross-sectional area occupied by the trees and is a stand density indicator commonly used by forest managers in silvicultural planning. The main objective of our study is to assess the benefit of combining ALS and optical satellite image time series (i.e., Sentinel-2) data for the downstream task of forest structure and biophysical variables estimation, in order to understand if multimodal information can improve the estimation of these physical quantities in a complex forest environment. While the former source of information, ALS, is widely used for the characterization of forest properties, the latter, Sentinel-2, brings knowledge about stand composition. Additionally, we also integrate in our analysis topographic information, which impacts (as confirmed by recent studies [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF] and empirically observed in the present experimental evaluation) the quality of both ALS structure information and Sentinel-2 radiometric signal.

To summarize, the main contributions of our research study are as follows:

1) the design of an end-to-end DL framework devoted to forest variables estimation via a multisource late-fusion approach; 2) a first study, to the best of our literature survey, that assesses the combination of topographic information, Li-DAR metrics, and Sentinel-2 multispectral and multitemporal data for total volume and basal area estimation in a complex forest environment. The rest of this article is organized as follows. Section II introduces the data available on the study site. Section III describes the proposed framework. The experimental settings and the results are reported and discussed in Sections IV and V, respectively. Finally, Section VI concludes this article.

II. DATA

A. Study Site

The study site is the Regional Natural Park of the "Massif des Bauges," located in France in the Alps mountains, between the two administrative departments of Savoie and Haute-Savoie (Fig. 1 and Fig. 2). The site is characterized by a steep and irregular topography, and an elevation ranging from 256 to 2217 m. The forest covers an area of 51 136 ha, which accounts for 60% of the total area of this inhabited park. The forest is, thus, a major component of the park's landscape and provides several ecosystem services to local inhabitants and neighboring cities: Wood production, outdoor activities, and biodiversity conservation. The forest stands consist of both deciduous and coniferous species, with a dominance of silver fir (Abies alba), Norway spruce (Picea abies), and common beech (Fagus sylvatica). Managed forests are mostly in an uneven-aged system.

B. Remote-Sensing Data

Airborne LiDAR data used in this study were not acquired under similar conditions over the whole area. Two surveys were conducted, each one covering a specific area. 1) The first acquisition (Acquisition 73) was conducted between June and September 2016, by the French National Institute of Geographic and Forest Information (IGN), and results from a collaboration between IGN and the governance and management board of Savoie Mont Blanc in order to provide the first LiDAR coverage over a whole French department in a mountainous area. On our study site, it covers an area of 53 600 ha in the southern part of the park, located in the Savoie department. Mean LiDAR point density is 4 points/m 2 . Raw data preprocessing, including geolocation, point clouds alignment, and data tiling was performed by a private contractor. 2) The second acquisition (Acquisition 74) was conducted in September 2018, by a private data provider in order to complete the coverage of the park area. It covers an area of 37 350 ha in the northern part of the park, located in the Haute-Savoie department. LiDAR point density is on average of 13 points/m 2 , with a strong spatial heterogeneity and local densities ranging from 4 to 25 points/m 2 over the whole area of acquisition. A 1-m spatial resolution digital terrain model (DTM) over the whole study site was derived from the LiDAR data. The LiDAR data of both acquisitions were normalized in height, using the DTM, to account for differences in terrain elevation.

Satellite data consists of a Sentinel-2 (S2) time series of five images, acquired between the years 2017 and 2018 at the following dates: October 7, 2017, November 1, 2017, August 28, 2018, September 27, 2018, and October 17, 2018. The images were chosen to ensure a cloud-free and snow-free coverage of the study site. Moreover, they were also selected to match, as much as possible, the dates of field and LiDAR surveys. The S2 images were obtained from the THEIA data platform 1 at level-2A (top of canopy reflectance). Ten spectral bands were considered in this study (Fig. 3) and all the bands were resampled (with the nearest neighbor strategy) to a 10-m spatial resolution. In addition, three topographic variables were computed at a spatial resolution of 10 m. The 1-m DTM was down-scaled to a 10-m raster of elevations, which was further used to compute both the aspect and slope layers, leading to a total of 53 channels describing the study site at a spatial resolution of 10 m.

C. Field Data

The GT data were obtained through a forest inventory campaign conducted by the French National Forest Office and the French National Research Institute for Agriculture, Food, and Environment between June and September 2018. A total of 291 circular plots, regularly distributed over the study site, were inventoried and georeferenced. Measurements were performed within concentric plots of 10-m and 15-m radius. Within the 15-m radius, all the trees with a diameter at breast height (DBH) above 17.5 cm had their DBH measured using a tape and the species were recorded. Within the 10-m radius, the trees between 7.5 and 17.5 cm DBH were classified and counted according to two criteria: 1) the diameter category (7.5 to 12.5 cm or 12.5 to 17.5 cm) and 2) the species category (coniferous or deciduous). Because no heights and not all the diameters were measured, the 1 [Online]. Available: https://www.theia-land.fr/en/homepage-en/

TABLE I ALS METRICS

database from the French National Forest Inventory produced by IGN was used to enrich the forest inventory with additional information to allow the computation of the two forest variables of interest in this study: Total volume and basal area. The variables were computed at the scale of a 15-m radius plot and successively converted to per hectare values. Fig. 4 shows the data distribution for each of the forest variables over the 291 plots.

D. Remote-Sensing-Based Variables

Inventory plots are formed by 30-m diameter disks. For each plot, the normalized LiDAR point cloud located within the extent of the 15-m radius circle was extracted. A threshold was applied to remove all the points below 5 m corresponding to lower vegetation. As shown in Table I, 55 standard LiDAR metrics were then computed at plot level using the R lidR package [START_REF] Roussel | LiDR: An R package for analysis of airborne laser scanning (ALS) data[END_REF]. Three topographic values were also averaged at plot level using the 1-m DTM, i.e., 1) elevation; 2) slope; and 3) aspect, leading to a total of 58 descriptors characterizing each plot. Regarding optical imagery, because of the spatial resolution of the resampled S2 images, the spectral information associated to each plot is extracted on a 3 × 3 pixel window that corresponds to a 30-m side square (referred to as a patch). The central pixel of a given patch is selected in order to contain the centre of the associated plot. The final 3 × 3 patch is then built around this central pixel. In addition to the 58 descriptors previously mentioned, each plot is also described by a 3 × 3 pixel patch of 53 channels at a 10-m spatial resolution extracted from the S2 time series plus features extracted from the 10-m DTM (see Subsection II-B).

To summarize, the labelled data used in this study consist of 291 plots described by two forest variables (total volume and basal area), 55 LiDAR-derived metrics, 3 topographic values at plot level, and 50 Sentinel-2 spectral bands plus 3 topographic bands at a 10-m spatial resolution. The details of all LiDAR-derived metrics and Sentinel-2 information are listed in Table I.

III. METHOD

In this section we describe our framework, named MMFVE, to cope with the fusion of multimodal data for the estimation of structure and biophysical variables. Firstly, we supply a general overview of our two-branch architecture and we detail the endto-end learning strategy we have adopted to learn the model parameters. Secondly, we introduce the per-source encoders we have designed to take into account the specificity of both LiDARderived metrics and Sentinel-2 time-series data.

A. Multimodal DL Architecture Fig. 5 depicts the proposed framework, MMFVE. In our scenario, each geospatial location is described by means of two different kinds of information: 1) metrics extracted from LiDAR acquisition and 2) spatial patches extracted from Sentinel-2 time series.

To manage and exploit the available multimodal information, MMFVE has two input branches, one for each of the sources to process. Each branch is associated to an encoder network that extracts a source specific representation: R S2 and R LIDAR . Subsequently, the different per-source representations are aggregated together considering a late-fusion schema [START_REF] Hong | More diverse means better: Multimodal deep learning meets remote-sensing imagery classification[END_REF] by concatenation of the per-source representations, with the aim to obtain a multimodal representation (R MM ), gathering together the multifaceted information describing the specific geospatial location, i.e., a 0.07 ha circular plot. Finally, the multimodal representation is fed through a hidden layer and a final output layer with only one neuron. The output neuron has no activation function associated to it and it only performs a linear combination of the input signals coming from the previous connected layer as commonly done in regression-based DL models [START_REF] Lathuilière | A comprehensive analysis of deep regression[END_REF]. The whole process is performed end-to-end.

The detailed description of the different components on which MMFVE is built on is reported in Fig. 5.

Regarding the general training process, since we are dealing with a regression task, we employ the Huber loss (HL) [START_REF] Lathuilière | A comprehensive analysis of deep regression[END_REF] as cost function. The HL is defined as follows:

HLoss(Y, Y ) = 1 2 (Y -Y ) 2 if |Y -Y | < Δ Δ|Y -Y | -1 2 Δ 2 otherwise (1)
where Y and Y are the true and the estimated values, respectively. Basically, the HL behaves as the mean squared error (MSE) loss when the error is lower than a threshold Δ, while it mimics the mean absolute error (MAE) loss otherwise. In addition, it is also differentiable at 0. The HL combines good properties from both MAE and MSE and avoids their limitations. Conversely to MAE, it bypasses large gradient backpropagation when the estimated quantity is getting closer to the real value and, differently from the MSE loss, it is more robust to outliers. In our scenario we set Δ to 1 for the HL.

B. Sentinel-2 and LiDAR Encoders

Regarding the S2 time-series data, we exploit a 2-D convolutional neural network (2-D-CNN) with the goal to take into account the spatial context describing the geospatial location associated to the forest variable to estimate. To this end, the final raster cube with 53 bands (50 Sentinel-2 bands and 3 bands containing topographic characteristics, as explained in Section II-D) is used as input. The encoder architecture has one 2-D convolutional layer with a kernel of size 3 × 3 and 256 filters, followed by two fully connected dense layers with 128 neurons each. The 2-D convolutional layer and the first fully-connected layers are followed by a dropout layer with a drop rate equal to 0.3.

Concerning the LiDAR information, describing the geospatial location of interest, we analyze the unstructured set of 55 metrics together with the additional three topographic variables (as described in Section II-D) by means of two fully-connected dense layers with 256 and 64 neurons. The first fully-connected layer is followed by a dropout layer with a drop rate equal to 0.3.

Once the encoder extracts per-source features (R S2 and R LiDAR ), they are concatenated (R MM ) and successively processed via a fully-connected dense layer with 256 neurons and ReLU activation function that, finally, feeds the output layer constituted by only one neuron.

IV. EXPERIMENTAL EVALUATION

In this section, we present the experimental settings and the results obtained on the study area introduced in Section II.

Regarding both forest variables (total volume and basal area), we firstly validate the importance of considering topographic information in addition to LiDAR-derived metrics for the estimation tasks. Then, we assess the use of only S2 and S2 plus topographic information. Successively, we evaluate the performance of our multimodal framework, MMFVE, to leverage the interplay between the different modalities we have as inputs (LiDAR, topographic, and S2 data) with respect to both machine and DL competitors. Finally, we report some information related to the training time associated to the different multimodal approaches.

As additional evaluations, we also assess the sensitivity of our framework to the loss function employed to optimize the internal network parameters (by leveraging standard loss functions for regression task [START_REF] Lathuilière | A comprehensive analysis of deep regression[END_REF]) and we supply a summary of the training time required by each of the employed competing methods, considering the full multisource setting (LiDAR metrics, topographic information, and S2 imagery).

A. Experimental Settings

Considering all the available input modalities in our multimodal setting, our first goal is to perform several preliminary experiments in order to better understand the contribution of each data source. We then evaluate the performances of our twobranch DL framework, MMFVE, with respect to its competitors.

Regarding LiDAR-derived metrics, as reported in Section II, each plot is described by 55 metrics plus three additional topographic variables (aspect, slope, and elevation). Concerning the S2 time series, enriched with the same three topographic variables, we consider input patches of size 3 × 3 (30 m × 30 m) that approximately correspond to the extent of a plot.

As competing approaches, we use standard and recent machine and DL techniques commonly employed to estimate forest variables at stand level from remote-sensing data [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF], [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF], [START_REF] Cosenza | Comparison of linear regression, k-nearest neighbour and random forest methods in airborne laser-scanning-based prediction of growing stock[END_REF], [START_REF] Coops | Modelling LiDAR-derived estimates of forest attributes over space and time: A review of approaches and future trends[END_REF]. More in detail, we adopt a K nearest neighbors regressor (KN N ) based on Euclidean distance, a random forest (RF ) based regressor method that leverages ensemble of trees via a bagging strategy to estimate the target variable, and the BART model recently evaluated in [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF]. In addition, we include as competing strateg,y the sparse stacked autoencoder (SSAE) model introduced in [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF]. This model is based on an early-fusion strategy to aggregate all the different information sources before the analysis step. It also provides an implicit comparison between early-and late-fusion strategies. According to [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF], the parameters of the SSAE model are learnt with MSE as loss function.

When only LiDAR-based (and topographic) information are considered, we include a multilayer perceptron (MLP ) with the same architecture as the encoder branch of MMFVE devoted to manage the LiDAR-derived metrics. Similarly, when only S2 imagery (and topographic) information are considered, we include a 2-D-CNN with the same architecture as the encoder branch of MMFVE devoted to manage the Sentinel-2 information.

The values of the different sources were normalized in the interval [0, 1]. The dataset was split into training, validation, and test set with 191, 50, and 50 plots, respectively. The evaluated models were optimized via training/validation procedure [START_REF] Ienco | Land cover classification via multitemporal spatial data by deep recurrent neural networks[END_REF]. For the KN N approach, the training/validation procedure allows to choose the best value of K in the set {1, 3, 5, 7, 9}, while for the RF regressor the optimization procedure permits to choose the best number of estimators (trees) in the set {200, 300, 400, 500, 600, 700, 800}, as well as the number of attributes for each subset of samples in the bagging procedure in the range [START_REF] Boyd | Satellite remote sensing of forest resources: Three decades of research development[END_REF][START_REF] Magnussen | Derivations of stand heights from airborne laser scanner data with canopy-based quantile estimators[END_REF]. For the BART method, we adopt the parameter setting employed in [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF].

Concerning MMFVE and the SSAE competitor, the learning stage was conducted over 400 epochs and the Adam optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] was used to learn trainable parameters with a learning rate of 10 -5 . Dropout rate was set to 0.3 and batch size was fixed to 1.

The assessment of the models' performances was done considering the test set and the following metrics: MAE, root mean squared error (rmse), and R 2 (R-squared, the coefficient of determination). The performance metrics are defined as follows:

MAE(Y, Y ) = 1 |Y | |Y | i=1 |Y i -Y i | (2) rmse(Y, Y ) = 1 |Y | |Y | i=1 (Y i -Y i ) 2 (3) R 2 (Y, Y ) = 1 - |Y | i=1 (Y i -Y i ) 2 |Y | i=1 (Y i -Ȳ ) 2 (4)
where Y contains the original values of the forest variable to estimate and Y is the set of values estimated by a particular approach.

Since model performances may vary depending on data splitting due to simpler or more complex samples involved in the different partitions, all metrics were averaged over 30 random splits of the dataset following the strategy mentioned above. The different neural network architectures were implemented using the Python Tensorflow library.

B. Results

Here, we provide the experimental evaluation according to different combinations of the input sources as well as the complete multimodal setting with LiDAR metrics, topographic information, and S2 imagery.

1) Assessing LiDAR Metrics With and Without Topographic Information: Tables II and III report the performances of the different competing approaches when LiDAR-derived metrics (L) and topographic information (T) are considered to estimate total volume and basal area, respectively. Source information Concerning the estimation of total volume (Table II), we can observe that the use of topographic information (in addition to LiDAR-derived metrics) generally ameliorates the results of the estimation algorithms. The only method that does not seem capable to exploit such additional information is the KN N algorithm. This can be explained by the fact that this approach relies on the computation of Euclidean distance and, to compute such a measure, all the input variables have the same importance. Since the number of topographic variables is much smaller than the number of LiDAR-derived features (3 versus 55), the contribution of the former is largely dominated by the contribution of the latter.

Conversely, RF , MLP , and BART are capable to leverage the useful information carried out by the topographic variables providing better performances when topographic information is considered (lower values of MAE and rmse and higher values of R 2 ). While RF and BART slightly improve their behavior, MLP clearly supplies the best performances gain regarding all the employed evaluation metrics. This latter method also provides the best absolute results outperforming all the other competing approaches when LiDAR-derived metrics alone and both LiDAR-derived and topographic variables are involved.

Regarding the estimation of basal area (Table III), we can observe a similar trend as the one depicted for the estimation of total volume. Also in this case, the KN N algorithm is unable of taking advantage of the extra information supplied by topographic variables, while the MLP approach exhibits the best performances with and without topographic variables (low MAE and rmse and high R 2 ), still demonstrating its ability to fully leverage the interplay between LiDAR-derived metrics and topographic information.

2) Assessing S2 Imagery With and Without Topographic Information: In order to understand the contribution of S2 time series in the estimation task, we have also assessed the performances of both forest variables (total volume and basal area) estimation considering only Sentinel-2 data with and without topographic information. The results are reported in Tables IV and V for the total volume and basal area forest variables, respectively. The achieved results are far inferior to The MLP (L, T ) method is also reported to assess the added value to also consider Sentinel-2 data.

the results previously reported, with the best results achieved when Sentinel-2 data are jointly exploited with topographic information via the 2-D-CNN model.

3) Multimodal Setting: Tables VI and VII summarize the performances of the different competing methods when all the multimodal information, describing the study area, are employed for the estimation of total volume and basal area, respectively. In this case, given a method (i.e., RF ), we indicate with {L,T,S2} the joint use of all the available information [i.e., RF (L, T, S2)]. In this comparison, as baseline to evaluate the added value of the Sentinel-2 time-series data, we consider the best performing method when only LiDAR-derived metrics and topographic information are used [MLP (L, T )].

For the case of total volume estimation (Table VI), we observe that MMFVE achieves the best results for all the evaluation metrics and the SSAE approach achieves the second best estimation. A direct comparison between our multimodal framework and MLP (L, T ) (that can be considered as an ablation of MM-FVE) underlines that the proposed neural network architecture is well suited to leverage the complementarity of the different input modalities. About the RF and the BART approaches, we can note that also in this case the use of Sentinel-2 data provides some improvement with respect to the results reported in Table II. But, we can also note that, in spite of feeding all the available data to the RF classifier, its performances are still lower than the performances achieved by the MLP (L, T ). Unlike the other methods, the use of all the available modalities negatively influenced the behavior of the KN N regression, thus, resulting in a systematic degradation of the evaluation metrics.

Concerning the estimation of the basal area variable, when all the available modalities are used as input, we observe similar behaviors as the ones exhibited for the estimation of total volume (Table VII). MMFVE outperforms all the other approaches demonstrating its ability to combine together multiple input modalities with SSAE, based on an early-fusion strategy, that achieves performances that are comparable to our framework. Both, the RF and the BART approaches, achieve marginal improvements with respect to its counterpart when Sentinel-2 input modality is missing. Table III highlights again their limited ability to fully exploit the interplay between LiDARderived metric, topographic information, and Sentinel-2 data. Also in this case, the KN N regressor is negatively affected by the large and heterogeneous amount of input information, emphasizing the inadequacy of this approach in a multimodal context.

Fig. 6 depicts the scatter plots of the estimated versus measured values for the two forest variables for two DL-based frameworks we have evaluated (MMFVE and MLP). More in detail, Fig. 6(a) and (b) represents the distributions related to total volume with only LiDAR and topographic information and with all the available sources (LiDAR, topographic, and Sentinel-2 data), respectively. We can observe that the latter configuration has higher estimation performances, and approximates the diagonal line, drawing the ideal scenario better than the former approach. This is also confirmed by the reported values for all the three evaluation metrics. Fig. 6(c) and (d) depicts the estimation of basal area when only LiDAR and topographic information are considered and when all the available sources (LiDAR, topographic, and Sentinel-2 data) are considered, respectively. Also in this case, we can see that the latter configuration, MMFVE, allows to achieve higher estimation performances. Similarly to what happens for the other forest variable, the use of all the available modalities permits to reach the best performances and results underline that the HL function allows some kind of improvement compared to both MSE and MAE regression loss functions. This is probably due to the fact that the HL combines good properties from both MAE and MSE and avoids their limitations: It bypasses large gradient backpropagation when the estimated quantity is getting closer to the real value, and it is more robust to outliers than MSE loss. 5) Time Performance Analysis for the Multimodal Setting: Finally, Table X summarizes the performances, in terms of training time, of all the different competing approaches in the case of multisource analysis (LiDAR metrics, topographic information, and Sentinel-2 imagery). We remind that all the experiments were performed on the Google Cloud Platform, Colab [START_REF] Bisong | Google Colaboratory[END_REF]. We can observe that the training times of all the methods are quite reasonable and affordable (spanning from a few seconds to 10 min.), with the highest running time exhibited by the BART approach, that takes around 9 min. to train a regression model.

V. DISCUSSION

To summarize, our research study proposes a novel multisource DL framework that adopts late-instead of earlyfusion [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF] strategy for the estimation of two forest variables with the aim to assess the combination of LiDAR metric and multispectral/multitemporal Sentinel-2 data. In addition, to the best of our literature survey, this is the first attempt to combine such sources of information via modern neural network approaches for the estimation of total volume and basal area variables in a complex forest environment. The proposed DL framework exhibits convincing performances in an operational scenario characterized by multimodal information as well as a realistic amount of available training samples.

Firstly, we have noted that the joint use of topographic information, LiDAR-derived metrics, and Sentinel-2 time series improved the estimation of both total volume and basal area with respect to the case in which only LiDAR metrics or only Sentinel-2 data are employed. For LiDAR-based total volume models, MAE, and rmse were reduced by -3.8% and -10.1%, for the RF and MLP models, respectively, and R 2 was increased by 5.6% and 10.3% for the same models when topographic data were introduced. Improvements were even slightly greater for basal area models. This improvement might be due to the fact that topographic information provides a kind of correction effect that compensates for the distortion of tree architecture induced by height normalization of the point cloud. This distortion increases with slope and can be considered as a source of signal artefacts. Trees growing in mountainous areas can also have curved trunks due to snow and wind. The spatial distribution of LiDAR 3-D point clouds and the derived metrics can thus change according to both slope and aspect for the same kind of stand type. Using topographic information is helpful to deal with such differences in point clouds. In the same spirit, we noticed that, despite the radiometric correction of topographic effects applied to Sentinel-2 images, differences in reflectances due to changes in exposition are still visible. Adding slope and aspect as input data in the model gives an opportunity to the neural network to better cope with these residual differences.

Secondly, Sentinel-2 data bring information mainly related to stand composition. Sentinel-2 data alone are not the best source of information to predict forest variables that are primarily driven by tree and stand structure, as also underlined in [START_REF] Wittke | Comparison of two-dimensional multitemporal Sentinel-2 data with three-dimensional remote sensing data sources for forest inventory parameter estimation over a boreal forest[END_REF]. For example, when comparing the best Sentinel-2 based models to the best LiDAR-based models, including topographic data in both cases, MAE and rmse increased by 38.7% and 31.7%, respectively, and R 2 decreased by 28.1% for total volume.

However, adding Sentinel-2 data to LiDAR information led to a significant improvement in the prediction with the proposed DL framework. R 2 increased by 13.4% and by 17.7% for total volume and basal area, respectively and both MAE and rmse were reduced by more than 20% for the two variables. These results underline the complementarity between this information, as also underlined in the recent study proposed in [START_REF] Ahmadi | Comparison of machine learning methods for mapping the stand characteristics of temperate forests using multi-spectral sentinel-2 data[END_REF], where the combination of the Sentinel-2 and POLSAR-derived topographic information improved the estimation of forest variables in temperate forests. At plot level, for a given volume or basal area, the 3-D distribution of LiDAR points and the derived metrics depend on tree species that are present in the plot. Therefore, the relationship between LiDAR-derived metrics and the targeted forest variables is dependant on stand composition. Results seem to demonstrate that this dependence can be addressed using Sentinel-2 data within an appropriate modeling framework. As a result, the quality of the predictions is unexpectedly high with regards to the complexity of the forest in the study area.

Thirdly, from a methodological point of view, when comparing competing approaches, we have observed that the MLP model exhibits the best behavior, in terms of evaluation metrics, for the analysis of both LiDAR-derived metrics and topographic information. We can explain this fact by the intrinsic ability of this approach to manage the different input information permitting to extract effective knowledge for the downstream regression task.

Fourthly, we have demonstrated the quality of our DL-based framework, based on late-fusion strategy, for the task of structure and biophysical variables estimation from heterogeneous (multimodal) remote-sensing data. The comparison with the strategy proposed in [START_REF] Zhang | Deep learning based retrieval of forest aboveground biomass from combined LiDAR and landsat 8 data[END_REF], based on an early-fusion option, advocates in favour of our framework, underlying the appropriateness of a multibranch architecture for our multimodal task. This result is completely in line with recent studies [START_REF] Hong | More diverse means better: Multimodal deep learning meets remote-sensing imagery classification[END_REF] that pinpoint the benefit of late-fusion approaches over early-fusion ones for multimodal remote-sensing analysis. More in detail, the obtained findings highlight that the availability of complementary remote-sensing data is not a direct synonym of performance improvements, but the way in which the diverse remote-sensing sources are combined together is crucial to fully exploit the interplay among the different input sources. This point is underlined by the empirical results we have reported in Section IV that pinpoint the ability of MMFVE to get the most out of the different information we can access.

Finally, we remind that our task is characterized by operational/realistic constraints, thus, resulting in a limited amount of GT information from which the relationship between multimodal remote-sensing data and the target forest variables is learnt. Despite such a data paucity setting, our approach clearly outperforms standard (i.e., KNN and RF regressor) and recent (i.e., BART and SSAE) competing strategies, demonstrating the quality of modern DL-based methods to provide competitive results also in scenarios where the amount of labelled samples is scarce. Dealing with data from multiple sources can also bring up additional issues due to the specificity of each available modality. Indeed, in this study, LiDAR and Sentinel-2 data were acquired through airborne and spaceborne configurations, respectively, which results in differences in terms of scale and resolution between these data sources. Moreover, the data georeferencing is of paramount importance in remote-sensing applications that rely on GT and reference data. In the context of multimodal data fusion, the input data on which models are trained are derived from multiple datasets, each acquired under specific conditions and preprocessed independently, and thus, misregistration among all the available modalities can negatively impact the final results. In this scenario, our proposed framework dealt with each modality independently through dedicated branches before performing a late-fusion step. This architecture allowed to take into account the peculiarities of each input data source and to process them in an adequate manner considering their differences in scale and resolution, before leveraging the interplay between all the sources in the internal-fusion step.

VI. CONCLUSION

In this work, we have evaluated the complementarity of Li-DAR, topographic, and Sentinel-2 information for the estimation of two forest variables (total volume and basal area) to characterize temperate forest properties in a multimodal scenario. To this end, we have designed and deployed a strategy named MMFVE, based on a DL framework, to cope with the richness and complexity of the available multimodal input data.

MMFVE is based on a two-branch architecture with dedicated per-source encoders and a latefusion step based on the concatenation of the per-source features. Thanks to an end-toend learning strategy, our framework is capable to effectively manage the available multimodal information associated to a particular spatial location.

The experimental evaluation with competing approaches, usually employed to cope with structural and biophysical variables estimation, highlights that the use of multimodal information alone does not guarantee an improvement of the final performances, but an effective exploitation of multimodal information tightly depends on the way the data-fusion process is conducted. More in detail, the results we have obtained on the estimation of total volume and basal area forest variables, in terms of all the considered evaluation metrics, underline the adequateness and the added value of our neural network based approach to fuse together LiDAR, topographic, and Sentinel-2 data.

Possible future works related to our research study can be devoted to a in-depth analysis of the interplay between LiDAR metrics and Sentinel-2 imagery on one hand, and topographic information on the other hand. For instance, the intensity channel derived from the laser wavelength of the LiDAR acquisition can be considered as an additional spectral band with respect to the original Sentinel-2 spectral bands due to the fact that it registers a different wavelength of the spectrum. Furthermore, evaluating which is the more appropriate way to integrate topographic information (at raster or metric level) as an additional and/or independent source of information (i.e., considering a threebranch architecture) can be another follow-up of the proposed research work. Finally, another possible avenue of research could be related to the analysis of the importance of the different input information in order to establish which LiDAR metric or spectral band the neural network retains most useful following a similar approach to the one proposed in [START_REF] Campos-Taberner | Understanding deep learning in land use classification based on Sentinel-2 time series[END_REF].

Fig. 1 .

 1 Fig. 1. Location of the study area in France and RGB composite from SENTINEL-2 imagery (top); Extent of the Massif des Bauges Natural Regional Park and location of inventory plots (Bottom).

Fig. 2 .

 2 Fig. 2. Location of the study area in France and RGB composite from SENTINEL-2 imagery (top); Extent of the Massif des Bauges Natural Regional Park and location of inventory plots (Bottom).

Fig. 3 .

 3 Fig. 3. S2 spectral bands used in this study are surrounded by a dashed green line.

Fig. 4 .

 4 Fig. 4. Distribution of GT data for the two forest variables basal area and total volume.

Fig. 5 .

 5 Fig. 5. Overview of MMFVE framework. The architecture has two branches, one devoted to the analysis of Sentinel-2 time-series data and one dedicated to the analysis of metrics extracted from LiDAR measurements. Then, the per-source feature representations are aggregated, under late-fusion strategy, by the means of the concatenation operation in order to perform the final forest-variable estimation.

Fig. 6 .

 6 Fig. 6. Estimated verus measured values for the total volume and basal area biophysical forest variables when only LiDAR and topographic information or all the available sources (LiDAR, topographic, and Sentinel-2 data) are considered. (a) Estimated versus measured total volume using LiDAR and topographic data by MLP (L, T ). (b) Estimated versus measured total volume using LiDAR, topographic, and S2 data by MMFVE. (c) Estimated versus measured basal area using LiDAR and topographic data by MLP (L, T ). (d) Estimated versus measured basal area using LiDAR, topographic, and S2 data by MMFVE.

TABLE II PERFORMANCES

 II OF THE DIFFERENT COMPETING APPROACHES (KNN, RF , BART , AND MLP ) IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN LIDAR-BASED INFORMATION IS CONSIDERED ALONE OR COUPLED WITH TOPOGRAPHIC INFORMATION TO ESTIMATE TOTAL VOLUME

TABLE III PERFORMANCES

 III OF THE DIFFERENT COMPETING APPROACHES (KNN, RF , BART , AND MLP ) IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN LIDAR-BASED INFORMATION IS CONSIDERED ALONE OR COUPLED WITH TOPOGRAPHIC INFORMATION TO ESTIMATE BASAL AREAis specified as input of the method. For instance, RF (L) indicates an RF taking as input only LiDAR-derived metrics, while RF (L, T ) indicates the same method with LiDAR-derived metrics and corresponding topographic variables as inputs.

TABLE IV PERFORMANCES

 IV OF THE COMPETING APPROACHES IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN SENTINEL-2 DATA ARE CONSIDERED ALONE OR COUPLED WITH TOPOGRAPHIC INFORMATION TO ESTIMATE TOTAL VOLUME TABLE V PERFORMANCES OF THE COMPETING APPROACHES IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN SENTINEL-2 DATA ARE CONSIDERED ALONE OR COUPLED WITH TOPOGRAPHIC INFORMATION TO ESTIMATE BASAL AREA TABLE VI PERFORMANCES OF THE DIFFERENT COMPETING APPROACHES IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN LIDAR-BASED INFORMATION IS JOINTLY EXPLOITED WITH SENTINEL-2 DATA AND TOPOGRAPHIC INFORMATION TO ESTIMATE TOTAL VOLUME

The MLP (L, T ) method is also reported to assess the added value to also consider Sentinel-2 data.

TABLE VII PERFORMANCES

 VII OF THE DIFFERENT COMPETING APPROACHES IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN LIDAR-BASED INFORMATION IS JOINTLY EXPLOITED WITH SENTINEL-2 DATA AND TOPOGRAPHIC INFORMATION TO ESTIMATE BASAL AREA

TABLE VIII PERFORMANCES

 VIII OF THE PROPOSED MMFVE FRAMEWORK IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN DIFFERENT LOSS FUNCTIONS ARE USED IN THE TRAINING PROCESS TO ESTIMATE TOTAL VOLUME

TABLE IX PERFORMANCES

 IX OF THE PROPOSED MMFVE FRAMEWORK IN TERMS OF MAE, RMSE, AND R 2 EVALUATION METRICS WHEN DIFFERENT LOSS FUNCTIONS ARE USED IN THE TRAINING PROCESS TO ESTIMATE BASAL AREA to make a step closer to the ideal scenario represented by the diagonal line portrayed in the scatter plot. 4) Evaluating Different Loss Functions for MMFVE: In addition, we compare the performances of our framework, MM-FVE, coupled with standard loss functions commonly employed for regression tasks: MSE and MAE. We indicate with MM-FVE(MSE), MMFVE(MAE), and MMFVE(HL) the MMFVE framework coupled with MSE, MAE, and HL functions, respectively. Tables VIII and IX depict the obtained results for total volume and basal area variables, respectively. The obtained

TABLE X RUNNING

 X TIMES FOR THE TRAINING PHASE OF MMFVE FRAMEWORK AND ITS COMPETITORS WHEN ALL THE AVAILABLE MULTIMODAL INFORMATION IS USED TO ESTIMATE ONE OF THE CONSIDERED FOREST VARIABLE

This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

This work was supported in part by the French Agency for Ecological Transition (ADEME) through the PROTEST project under Grant 1703C0069, in part by GRAINE program, and in part by the French Region of Occitanie with the funding of part of K. Dayal's Ph.D. scholarship.


S Leblanc

M Lefort

E Le Page

L Michel

Emmanuelle Leray
email: emmanuelle.leray@ehesp.fr

Trends in disease-modifying therapy use in patients with multiple sclerosis using a 10-year population-based cohort study in France

Keywords: disease-modifying therapies, health administrative database, multiple sclerosis, time-series, untreated

www.indsante.fr/).

Introduction

Multiple sclerosis (MS) is a neurological chronic disease starting in young adulthood. There is no cure but several disease-modifying therapies (DMTs) are used to reduce occurrence of relapses and delay disability progression. In recent years, a number of novel and increasingly effective DMTs have been approved [START_REF] Mcginley | Diagnosis and Treatment of Multiple Sclerosis: A Review[END_REF](. The availability of new treatment options for MS provides an opportunity for improving tolerance, drug adherence, quality of life and long-term clinical outcomes but it also makes patient management more complex. When choosing a drug, neurologists should take into consideration several parameters including patient and disease characteristics, mechanisms of action, effectiveness, safety profile, way of administration, pregnancy desire and patient preferences. For instance, first-generation immunomodulators (beta interferon and glatiramer acetate -injectable drugs) have well-established short-term and long-term safety profiles, while new oral first-line drugs such as teriflunomide and dimethyl fumarate may appear more convenient from the patient perspective.

Our team recently published a detailed description of the use of DMTs for MS at French national level using the exhaustive MS population issued from the national health data system over 2010-2015 (98% coverage) [START_REF] Leblanc | Disease-modifying therapy usage in patients with multiple sclerosis in France: A 6-year population-based study[END_REF]. First-line DMTs were the mostly used, particularly injectable ones. However, we do expect that therapeutic practices have changed since 2015, as highlighted by another French study using the same database about 10,240 patients initiating a first-line DMT between 2014 and 2016 [START_REF] Vermersch | An analysis of first-line disease-modifying therapies in patients with relapsing-remitting multiple sclerosis using the French nationwide health claims database from 2014-2017[END_REF], where oral DMTs had the predominant position (teriflunomide 34.6% and dimethyl fumarate 33.7%).

Moreover, based on our previous work on insurance claims data [START_REF] Leblanc | Disease-modifying therapy usage in patients with multiple sclerosis in France: A 6-year population-based study[END_REF], half of patients did not receive any DMT over the 6-year study period while it was estimated to be one third in 2016-2017 in MS expert centers [START_REF] Moisset | Untreated patients with multiple sclerosis: A study of French expert centers[END_REF]. Such discrepancy can be linked to either data source (i.e. recruitment) or study period (i.e. changes over time) or both.

In that context, we decided to use a large nationwide sample based upon 10-year administrative data to describe time trends in the utilization of DMTs in MS patients in France. The main objective was to describe changes in the use of DMTs over the period 2009-2018, by comparing initial and subsequent DMTs between four successive periods. The secondary objective was to study the impact of the arrival of oral DMTs on the use of injectable DMTs, by simulating what would have been the prescriptions without these new drugs.

Materials and methods

Data source

Data were extracted from the General Sample of Beneficiaries ('Echantillon Généraliste des Bénéficiaires', EGB). The EGB database is a representative 1/97 sample of the population (around 600,000 people) covered by the French national healthcare insurance system [START_REF] Tuppin | French national health insurance information system and the permanent beneficiaries sample[END_REF][START_REF] Martin-Latry | Pharmacoepidemiological research using French reimbursement databases: yes we can![END_REF][START_REF] Moulis | French health insurance databases: What interest for medical research?[END_REF][START_REF] Palmaro | Overview of drug data within French health insurance databases and implications for pharmacoepidemiological studies[END_REF]. It contains anonymized demographic characteristics and exhaustive records of healthcare reimbursements including drugs dispensed in retail pharmacies, as well as data on hospital admissions with their related International Classification of Diseases, 10th version (ICD-10) diagnoses.

Study population

The study population included all patients with MS identified in the EGB database over the study period, i.e. between 1 January 2009 and 31 December 2018. Patients were considered as having MS if they met at least one of the following criteria: (i) having at least one reimbursement for an MS-specific DMT (beta-interferon, glatiramer acetate, fingolimod, natalizumab, teriflunomide or dimethyl fumarate); (ii) having a long-term disease (LTD) registration for MS, or (iii) having at least one hospital admission with a diagnosis of MS based on the ICD-10 G35 diagnosis code [START_REF] Foulon | Prevalence and mortality of patients with multiple sclerosis in France in 2012: a study based on French health insurance data[END_REF]. In France, being registered with LTD allows individuals with chronic diseases to receive full healthcare reimbursement. In total, 30 diseases are currently included in the list of LTD, including MS.

The date of MS identification was defined as the earliest date between the date of LTD registration for MS (can be before 2009), the first date of hospitalization for MS or the first date of MS-specific DMT prescription from 2009 to 2018. Based upon data availability, study entry was set at respectively 1 January 2009 for prevalent cases (i.e. MS identified through LTD registration before 2009) and date of MS identification for incident cases (2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The follow-up period ended with death or the study end on 31 December 2018, whichever came first.

Data

Patients with MS were defined as treated if they had at least one reimbursement of the following drugs: beta-interferon, glatiramer acetate, fingolimod, natalizumab, teriflunomide, dimethyl fumarate, azathioprine, methotrexate, mycophenolate mofetil or rituximab. These last 4 drugs were not considered to identify MS cases as they are not MS-specific. Neither mitoxantrone, nor cyclophosphamide and alemtuzumab were considered since they were not available in the hospital data. Neither ocrelizumab nor cladribine or ozanimod were available in 2018 in France.

Demographic patient characteristics included gender, age at MS identification, LTD status (for MS and for other conditions), follow-up duration from MS identification, and occurrence of death. As cause of death was not available in the dataset, for in-hospital deaths, the diagnosis of the hospital stay was used to approximate cause of death.

Molecules and start/stop dates from drugs claims were considered to describe therapeutic management of MS. Initial DMT and subsequent ones were considered. DMTs were also classified into four groups: first-line injectable (beta-interferon and glatiramer acetate), first-line oral (dimethylfumarate and teriflunomide), second-line (natalizumab and fingolimod), and off-label (azathioprine, methotrexate, mycophenolate mofetil and rituximab). Untreated periods were also considered in the analyses, and defined as a minimum of 3 months without any drug claim.

Four periods of MS identification were defined (before 2009, 2009 to 2011, 2012 to 2015 and 2016 to 2018), in order to assess changes in MS therapeutic management. For prevalent cases (identified before 2009), age at study entry and age at MS identification were computed, while for incident cases (from 2009 to 2018), study entry always corresponded to MS identification.

Statistical analysis

Characteristics of patients and DMTs use were described from study entry, overall and for each time period. The time from study entry to death was calculated using Kaplan-Meier estimates [START_REF] Kaplan | Nonparametric estimation from incomplete observations[END_REF].

Characteristics of patients who took at least one DMT and characteristics of patients who did not receive any DMT over the study period were also described overall and for each time period. Prescriptions of each DMT were computed for each patient from study entry on a monthly basis and summarized in a chronogram. Such graph displays the overall pattern of therapeutic use, with the cumulative proportion of patients per DMT (y-axis) for each time unit (x-axis). In addition, therapeutic sequences, i.e. initial drug and subsequent ones, were represented using sunburst diagrams, as such graphs make possible to visualize hierarchical data, depicted by concentric circles. The circle in the centre represents the root node (the initial DMT in the present case), with the hierarchy moving outward from the center (subsequent DMTs). Regarding comparisons according to time periods, as theoretical follow-up duration decreased over time, several indicators were calculated on a duration truncated to two years, in order to make them comparable.

All analyses were performed using the SAS software version 9.4 (SAS Institute Inc., Cary, NC, USA) except time series analyses which were performed using the fpp2 library [12] of the R version 3.6.1 software [START_REF] Team | A language and environment for statistical computing[END_REF].

Results

Characteristics of the study population

In the EGB sample which represents 1/97 of the whole French health insurance database, 1,508 patients were identified as having MS over the period 2009-2018. Almost half of the patients (48.9%, n=738) were identified before 2009 (prevalence rate: 129 per 100,000) as they presented a LTD status for MS, and the remaining were incident cases over the study period. Indeed, 15.6% (n=236) were identified between 2009 and 2011, 20.9% (n=315) between 2012 and 2015 and 14.5% (n=219) between 2016 and 2018. Characteristics of the overall study population as well as subgroups according to periods of MS identification are presented in Table 1A. Median age at MS identification was 41 years and 70.2% of patients were women. Among prevalent cases, age at study entry was 49 years while age at MS identification was 39 years. Overall, one third of patients with MS had another long-term disease status, reflecting the burden of comorbidities. In addition, 10.7% of patients died over the median follow-up duration of 9 years. Among the 161 deaths, median age at death was 69 years, and 29.8% occurred at hospital; the medical diagnosis of this stay was MS in 77.1% of them. Moreover, 50.9% of the deceased patients (n=82) had at least another LTD status, cardiovascular diseases (n=58) and cancer (n=33) being the most common.

Comparison of treated with untreated patients

Overall, 876 patients (58.1%) were treated at least once over the study period. As shown in Table 1B, treated patients were younger than untreated patients. As expected, LTD registration for MS was more frequent among treated patients than untreated patients while LTD registration for another disease was higher among untreated patients, probably reflecting older age and comorbidities. In accordance, the crude survival probability was higher in patients who received at least one DMT over the study period than in untreated ones (Figure 1).

Use of DMTs over time

Use of DMTs over time is described on Figure 2 where monthly proportions of patients per DMT are presented. About two-thirds of patients were not treated each month and this proportion was stable over time. These repeated cross-sectional proportions have to be confronted with the proportion of patients who were never treated on the 10-year period, i.e. 41.9% (Table 1B).

The use of the new oral treatments as soon as they arrived was clear, with fingolimod from 2011 (orange) and dimethyl fumarate and teriflunomide from 2014 (light-blue). Marketing of the two first-line oral drugs has led to a large decrease in the use of first-line injectable drugs, especially beta-interferon (darkblue), while second-line and off-label drugs did not seem to be impacted. At the end of the study period, it should be noted that all products seemed to share almost equal parts (beta interferon, teriflunomide, dimethyl fumarate, and fingolimod, then natalizumab and glatiramer acetate). This shift to new marketed drugs was confirmed using the time-series model, as shown in Figure 3, when comparing the expected prescriptions (dotted lines) to the actual ones (solid lines) on the period 2014-2018.

Therapeutic sequences over time

Additional characteristics of longitudinal therapeutic sequences are summarized in Table 2 and Figure 4. Initial DMT was mostly first-line DMTs (around 75-80%), followed by off-label (around 10-15%) and second-line DMTs (around 10%). Nearly half of patients switched from initial DMT to another DMT over the follow-up, and when restricting follow-up duration to the first two years, the proportion of patients receiving more than one DMT seemed to increase from ancient periods to more recent periods. In 2012-2015, the proportion rose to 27.6%, probably linked to the arrival of oral first-line drugs. Overall, 24.3% of patients received a second-line DMT, without period effect. Sunburst diagrams (Figure 4) showed the variety of therapeutic sequences, overall and in each time period. Indeed, many patients received different kinds of DMTs, as represented by the different colors in the concentric bands. Here again, the shift from injectable to oral first-line drugs was clear over time.

By contrast, use of second-line DMTs (in red) either in escalating strategy or as initial treatment seemed stable over time. We can notice that patients also regularly underwent discontinuation of DMTs (in yellow).

Discussion

The present study aimed to describe temporal trends in therapeutic management of MS in France over the period 2009-2018. As expected, we showed that patients were treated earlier over time and many of them received several DMTs. Moreover, in the past 10-year period, it seems that as soon as a drug becomes available, it is used to the detriment of previous ones, leading to a strong shift from injectable to oral first-line DMTs over time.

The data were based on a representative 1/97 random sample of the French population covered by the health insurance system (98%), which reduces selection bias and offers the opportunity to describe the wide MS population. The characteristics (e.g. gender, age, and LTD status) of the study population are very close to those of MS patients in France in general [START_REF] Roux | Use of healthcare services by patients with multiple sclerosis in France over 2010-2015: a nationwide population-based study using[END_REF].

Overall, our results are in line with a recent Canadian descriptive study [15] where people who were not treated during follow-up were older and had a higher comorbidity burden than people who were treated.

Nonetheless, it can be noticed that nearly 1 in 4 treated patients had at least some comorbidity, and/or were 50 years or older.

Our results regarding DMTs use are also consistent with those previously published about DMT initiations in 2016-2018 [START_REF] Vermersch | An analysis of first-line disease-modifying therapies in patients with relapsing-remitting multiple sclerosis using the French nationwide health claims database from 2014-2017[END_REF]. Indeed, we found that the proportion of patients treated with beta interferon (initiations and follow-up) decreased by half when new oral drugs became available. At the end of the study period, relative parts of each first-line DMT were almost similar. Our results are also in accordance with a French nation-wide study based upon the insurance regimen dedicated to self-employed workers (6% of the French general population) which found an overall DMT utilization of 34% in 2013 and 38% in 2015 [16]. They also showed that the first-generation DMT utilization rate decreased to the profit of second-generation drugs.

We also showed that about two-thirds of patients with MS were not treated each month, and 41.9% were never treated during the 10-year period. It should be noted that never treated patients in the present study do not necessarily correspond to patients who never received DMT throughout the whole disease course, as the term "never" only refers to the limited observation period, i.e. 2009-2018. Such estimates are consistent with annual untreated proportions of 66% and 61% in 2013 and 2015, respectively [16], but are higher than the 31% of patients from MS expert centers who were not receiving a DMT in 2016-2017 [START_REF] Moisset | Untreated patients with multiple sclerosis: A study of French expert centers[END_REF]. When looking at international literature, 29% of MS cases filled a DMT prescription in Canada during the 22-year study period [15], 10% to 50% in Denmark depending on MS onset period [17], half of patients enrolled in a private health insurance program versus one third of those enrolled in a public program in the US [18] and 80% in MSBase [19]. The differences in the proportions of untreated between studies likely reflect differences in patients' selection (age, MS phenotype, health status…), study period, and healthcare systems. Indeed, use of DMT was generally higher in the more recent periods in line with the increased availability of drugs, as well as in patients selected from MS expert centers, such as OFSEP [START_REF] Moisset | Untreated patients with multiple sclerosis: A study of French expert centers[END_REF] or MSBase [19] series. In population-based studies ([15;17]; the present one), there is no recruitment bias, which may represent a significant advantage over recruitment through MS centers which focus on patients who are actively seeking expert care and may led to more severe MS patients [20]. In France, DMTs are fully reimbursed by the health insurance system, which is not the case in many countries. We could assume that the absence of personal expense (out-of-pocket cost) should facilitate access to treatment, but obviously, it does not mean that the proportion of untreated patients is decreased. As the untreated subgroup was found to be older than the treated one, it probably included patients with progressive MS, but information related to MS phenotype was not available in the present dataset.

Among treated patients, we also found regular periods without treatments which may reflect different situations: treatment discontinuation due to side effects, treatment discontinuation before and/or during pregnancy, washout period required when switching from one agent to another, other personal reasons (convenience, lassitude, travel…). This was also highlighted in patients coming from MS expert centers, where treatment exposure was only 60% of total follow-up [START_REF] Moisset | Untreated patients with multiple sclerosis: A study of French expert centers[END_REF].

The variety in therapeutic sequences is in line with findings previously published by our team [START_REF] Leblanc | Disease-modifying therapy usage in patients with multiple sclerosis in France: A 6-year population-based study[END_REF] and illustrates both the lack of clear and specific recommendations in France and the shared indications among drugs. It also reflects the increased number of available drugs, as well as the expanding literature about the comparative effectiveness and safety in real-word settings, and therapeutic strategies. Overall, half of patients were exposed to only one DMT over follow-up, nearly 80% of people filled a prescription for a first-line DMT and one quarter for a second-line DMT.

In the last 20 years, the treatment options for relapsing-onset MS have considerably been enriched with a variety of new drugs in order to reach a better control of the disease. In that context, neurologists are asked to initiate effective treatment early in the disease course, monitor treatment response, and modify or switch therapies in case of failure. We can assume that the level of expectations regarding disease control is increasing from both neurologists and patients, and that a new goal of "no evidence of disease activity" (NEDA) may be supplanting the previous aim of relapse rate reduction [21]. However, in the present work, we were not able to measure effectiveness of the therapeutic strategies, as our main limitation was the lack of clinical data. Indeed, we observed therapeutic practices from French neurologists and changes over time, but we cannot assess whether the DMT choices were appropriate, whether switch decisions were done at the right time, or whether untreated patients did really not need to be treated, as we did not have information regarding relapses, disability, MRI parameters. Even MS phenotype was not available, meaning that we could not stratify findings between relapsing and progressive forms. In addition, the last period of MS identification (2016-2018) was probably too short to reflect current practices. However, such descriptive real-world studies can offer the opportunity to describe use of DMTs in subpopulations that are generally excluded from randomised clinical trials, such as old patients or late-onset MS or pregnant women.

Accepted manuscript / Final version

Accepted manuscript / Final version

Accepted manuscript / Final version

  Then, time series of each DMT from January 2009 until February 2014 (first first-line oral drug availability) were used to forecast trends from March 2014 to December 2018, to simulate what would have happened in the absence of new marketed oral drugs. For that purpose, a seasonal autoregressive integrated moving average (SARIMA) model [11] was implemented. The component models were fitted independently on the training period and then used to forecast the expected numbers over the test period. For each component model, 80% prediction intervals were estimated based on the standard deviation of each step forecast [12].

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figures Legend

Fig1 Accepted manuscript / Final version 20 Fig2

 20 Fig1

Table 1A .

 1A Characteristics of the study population, overall and according to periods of multiple sclerosis (MS) identification (n=1,508)

		Overall cohort	Before 2009	2009-2011	2012-2015	2016-2018
	Period of MS identification
		N=1,508	n=738 (48.9%)	n=236 (15.6%)	n=315 (20.9%)	n=219 (14.5%)
	Women, n(%)	1,058 (70.2%)	540 (73.2%)	155 (65.7%)	218 (69.2%)	145 (66.2%)
	Age at MS identification (years)*	41 (32-50)	39 (31-48)	43 (34-54)	40 (30-51)	43 (32-53)
	Age at study entry (years)*	45 (36-57)	49 (40-59)	43 (34-54)	40 (30-51)	43 (32-53)
	Time from MS identification in 2009 (years)*	0 (0-8)	8 (4-13)	--	--	--
	Long-term disease for MS at entry, n(%)	1.278 (84.7%)	738 (100.0%) 162 (68.6%)	222 (70.5%)	157 (71.7%)
	Age at LTD for MS registration (years)*	40 (31-48)	39 (31-48)	42 (33-49)	39 (30-48)	42 (31-49)
	At least another LTD over the study period, n(%)	509 (33.8%)	229 (31.0%)	97 (41.1%)	111 (35.2%)	72 (32.9%)
	Follow-up time (years)*	9 (4-10)	10 (10-10)	9 (8-9)	5 (4-6)	2 (1-2)
	Deaths over the study period, n(%)	161 (10.7%)	104 (14.1%)	33 (14.0%)	21 (6.7%)	3 (1.4%)

*median (q1-q3); MS: Multiple Sclerosis, LTD: Long-Term Disease.

Table 1B .

 1B Characteristics of patients treated at least once and untreated patients over the study period, overall and according to periods of multiple sclerosis (MS) identification (n=1,508)

	Period of MS	Overall cohort	Before 2009	2009-2011	2012-2015	2016-2018
	identification	N=1,508	n=738 (48.9%)	n=236 (15.6%)	n=315 (20.9%)	n=219 (14.5%)
		Treated	Not	Treated	Not	Treated	Not	Treated	Not	Treated	Not
			treated		treated		treated		treated		treated
		n=876	n= 632	n=429	n=309	n=134	n=102	n=196	n=119	n=117	n=102
		(58.1%)	(41.9%)	(58.1%)	(41.9%)	(56.8%)	(43.2%)	(62.2%)	(37.8%)	(53.4%)	(46.6%)
	Women, n(%)	619 (70.7%)	439 (69.5%)	319 (74.4%)	221 (71.5%)	88 (65.7%)	67 (65.7%)	133 (67.9%)	85 (71.4%)	79 (67.5%)	66 (64.7%)
	Age at MS
	identification
	(years)*
	LTD for MS at entry,	817 (93.3%)	461 (72.9%)	429 (100.0%)	308 (99.7%)	118 (88.1%)	44 (43.1%)	171 (87.2%)	51 (42.9%)	99 (84.6%)	58 (56.9%)
	n(%)
	At least another	232 (26.5%)	277 (43.8%)	111 (25.9%)	118 (38.2%)	43 (32.1%)	54 (52.9%)	53 (27.0%)	58 (48.7%)	25 (21.4%)	47 (46.1%)
	LTD, n(%)
	Follow-up time	10 (5-10)	8 (3-10)	10 (10-10)	10 (10-10)	9 (8-10)	8 (7-9)	5 (4-6)	4 (4-6)	2 (1-3)	1 (1-2)
	(years)*
	Deaths over the	40 (4.6%)	121 (19.1%)	27 (6.3%)	77 (24.9%)	8 (6.0%)	25 (24.5%)	5 (2.6%)	16 (13.4%)	0 (0.0%)	3 (2.9%)
	study period, n(%)

*median (q1-q3); MS: Multiple Sclerosis, LTD: Long-Term Disease.

Table 2 .

 2 Characteristics of the use of disease modifying therapies (DMTs) in patients with multiple sclerosis who were treated at least once over the study period (n=876)

	Overall cohort	Before 2009	2009-2011	2012-2015	2016-2018
	n= 876 (100.0%)	n=429 (49.0%)	n=134 (15.3%)	n=196 (22.4%)	n=117 (13.4%)

*median (q1-q3); MS: Multiple Sclerosis; DMT: Disease Modifying Therapy.

the authors cannot publicly release the data from the French national health data system (SNDS). However, a request for data reuse may be made and would require the obtention of prior approval of the French regulatory authorities (https://www.snds.gouv.fr/SNDS/Processus-d-acces-aux-donneesandhttps://