domain: sdv stat math

text:
Yoonjee Kang

Denis Thieffry

Laura Cantini
email: cantini@bio.ens.psl.eu

Evaluating the reproducibility of single-cell gene regulatory network inference algorithms

Keywords: biological networks, scRNA-seq, single-cell, transcriptome, network inference

Networks are powerful tools to represent and investigate biological systems. The development of algorithms inferring regulatory interactions from functional genomics data has been an active area of research. With the advent of single-cell RNA-seq data (scRNA-seq), numerous methods specifically designed to take advantage of single-cell datasets have been proposed. However, published benchmarks on single-cell network inference are mostly based on simulated data. Once applied to real data, these benchmarks take into account only a small set of genes and only compare the inferred networks with an imposed ground-truth. Here, we benchmark six single-cell network inference methods based on their reproducibility, i.e. their ability to infer similar networks when applied to two independent datasets for the same biological condition. We tested each of these methods on real data from three biological conditions: human retina, T-cells in colorectal cancer, and human hematopoiesis. Once taking into account networks with 100,000 links, GENIE3 results to be the most reproducible algorithm and, together with GRNBoost2, show higher intersection with ground-truth biological interactions. These results are independent from the single-cell sequencing platform, the cell type annotation system and the number of cells constituting the dataset. Finally, GRNBoost2 and CLR show more reproducible performance once a more stringent thresholding is applied to the networks (1000-100 links). In order to ensure the reproducibility and ease extensions of this benchmark study, we implemented all the analyses in scNET, a Jupyter notebook available at https://github.com/ComputationalSystemsBiology/scNET.

Introduction

Biological systems are inherently complex, in particular because of the emergent phenotypic properties arising from the interaction of their numerous molecular components. Characterizing genotype to phenotype connections and pathological deregulations thus requires to identify the biological macromolecules involved (e.g. genes, mRNAs, proteins), but also how these interact in a huge diversity of cellular pathways and networks [START_REF] Barabási | Network biology: understanding the cell's functional organization[END_REF]).

In the post-genomic era, biological networks have been extensively exploited to investigate such complex interactions among biological macromolecules [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF]Sonawane et al., 2019;[START_REF] Silverman | Molecular networks in Network Medicine: Development and applications[END_REF]. Network-based studies brought crucial insights into cell functioning and This is a provisional file, not the final typeset article

Running Title

diseases [START_REF] Basso | Reverse engineering of regulatory networks in human B cells[END_REF][START_REF] Margolin | ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF][START_REF] Ideker | Protein networks in disease[END_REF]. A network is a graphbased representation of a biological system, where the nodes represent objects of interest (e.g. genes, mRNAs, proteins), while the edges represent relations between these objects (e.g. gene co-expression, or binding between two proteins). Different approaches can be used to reconstruct biological networks.

Here, we focus on data-driven methods, which infer networks from gene expression data with the help of reverse engineering techniques (Sonawane et al., 2019). Network inference algorithms were first proposed to extract information from bulk gene expression data, and their development has been an active area of research for more than 20 years [START_REF] Barabási | Network medicine: a network-based approach to human disease[END_REF]The DREAM5 Consortium et al., 2012;Verny et al., 2017;Sonawane et al., 2019;[START_REF] Silverman | Molecular networks in Network Medicine: Development and applications[END_REF]. With the advent of single-cell RNA sequencing (scRNA-seq), we started to gather transcriptomic data from individual cells, enabling proper studies of their heterogeneity. However, the analysis of scRNA-seq data comes with a variety of computational challenges (e.g. small number of sequencing reads, systematic noise due to the stochasticity of gene expression at single-cell level, dropouts) that distinguish this data type from its bulk counterpart. For this reason, network inference methods originally developed for bulk gene expression data may not be suitable for data generated from single cells. The development of network inference algorithms has thus recently undergone a strong shift towards the design of methods targeting single-cell data [START_REF] Fiers | Mapping gene regulatory networks from single-cell omics data[END_REF].

Two benchmarks of single-cell network inference methods have been published [START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF][START_REF] Pratapa | Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[END_REF]. Both works evaluate network inference algorithms by comparing the inferred network with a ground-truth. These works are also mostly focused on simulated data and they apply a strong filtering on genes (leaving only 100-1,000 genes for network inference). Chen et al. [START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF] considered five methods targeting bulk data and three methods specifically designed for single-cell data. More recently, Paratapa et al. [START_REF] Pratapa | Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[END_REF] focused on twelve methods designed for single-cell data. Both benchmarks concluded that the overall performances of all methods were quite disappointing, and that network inference remains a challenging problem.

Here, we evaluate network inference algorithms based on their reproducibility, i.e. their ability to infer similar networks once applied to two independent datasets for the same biological condition (e.g. two independent scRNA-seq datasets obtained from colorectal tumors). The rationale behind this comparison is that, if the two independent datasets are profiled from the same biological condition (e.g. colorectal cancer) involving the same cell types, we can expect that the regulatory programs underlying them should strongly overlap. As a consequence, a good network inference algorithm should infer highly overlapping networks when applied to single-cell datasets profiled from the same biological condition. We selected six algorithms spanning the main network inference formulations that do not require an ordering of the cells according to pseudo-time, and we tested the reproducibility of the inferred networks in three biological systems: human retina, T-cells in colorectal cancer and human hematopoiesis. Differently from previous benchmarks, we only applied a soft filtering on genes, thus testing the algorithms based on their performances to infer networks involving from 6,000 to 12,000 nodes/genes. From our benchmark, once an high number of links is taken into account (100,000), GENIE3 results to be the most reproducible algorithm and, together with GRNBoost2, show the highest intersection with ground-truth biological interactions. GRNBoost2 and CLR have instead better performances for low link numbers (1000-100). In order to ensure the reproducibility and ease extensions of this benchmark study, we implemented all the analyses in a Jupyter notebook, called scNET and available at https://github.com/ComputationalSystemsBiology/scNET . Materials and Methods

2.1

Benchmarked single-cell network inference algorithms Starting from the exhaustive collection of single-cell network inference algorithms presented in [START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF][START_REF] Pratapa | Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[END_REF], two main categories of methods can be distinguished. Some methods interpret scRNA-seq as time-course expression data, where the pseudo-time corresponds to the time information. These methods are frequently based on Ordinary Differential Equations (ODEs) and are relevant for biological systems undergoing dynamic transcriptional changes (e.g. scRNA-Seq performed on differentiating cells) [START_REF] Matsumoto | SCODE: an efficient regulatory network inference algorithm from single-cell RNA-Seq during differentiation[END_REF]. In contrast, other methods do not use pseudo-time information to infer networks. These methods generally use statistical measures (e.g. correlation, mutual information) to infer regulatory connections and are thus better suited for transcriptomic data not affected by strong dynamical processes (e.g. retina cells in normal state). Testing reproducibility strictly requires the availability of two independent scRNA-seq datasets reflecting the same biological condition and presenting as few as possible technical variations. Indeed, the presence of technical variations due to the sequencing or experimental procedures could drastically impact the outcome of our comparison. In this respect, finding independent scRNA-seq datasets reflecting dynamic transcriptional changes, generated with the same experimental procedure, is really challenging. We thus decided to focus our benchmark study on network inference methods that do not use the pseudo-time information. In addition, only algorithms provided in R or Python code are here taken into account. Six single-cell network inference methods are thus considered in this evaluation: GENIE3 [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF], GRNBoost2 [START_REF] Moerman | GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks[END_REF], PPCOR [START_REF] Kim | ppcor: An R Package for a Fast Calculation to Semi-partial Correlation Coefficients[END_REF], PIDC [START_REF] Chan | Gene Regulatory Network Inference from Single-Cell Data Using Multivariate Information Measures[END_REF], CLR [START_REF] Faith | Large-Scale Mapping and Validation of Escherichia coli Transcriptional Regulation from a Compendium of Expression Profiles[END_REF] and GeneNet [START_REF] Opgen-Rhein | From correlation to causation networks: a simple approximate learning algorithm and its application to high-dimensional plant gene expression data[END_REF]. All the methods selected for this benchmark were originally designed for bulk data and they span the main mathematical formulations of network inference, as described in (The DREAM5 Consortium et al., 2012). Of note, GENIE3, GRNBoost2 and PIDC are also the best performing in the single-cell benchmark of Pratapa et al. GEne Network Inference with Ensemble of Trees (GENIE3) [START_REF] Huynh-Thu | Inferring regulatory networks from expression data using tree-based methods[END_REF]) is a tree-based network inference method. For each gene G1 in the expression dataset, GENIE3 solves a regression problem, determining the subset of genes whose expression is the most predictive of the expression of G1. This method was the best performing algorithm in the DREAM4 In Silico Multifactorial challenge [START_REF] Greenfield | DREAM4: Combining genetic and dynamic information to identify biological networks and dynamical models[END_REF]. GENIE3 requires in input the scRNA-seq expression matrix and a list of Transcription Factors (TFs). In our tests the list of human TFs provided in input corresponds to the intersection between the expressed genes and those annotated as encoding TFs by [START_REF] Chawla | TFcheckpoint: a curated compendium of specific DNA-binding RNA polymerase II transcription factors[END_REF]. The output of GENIE3 is a weighted network linking TFs with predicted target genes. The weight associated with each link corresponds to its Importance Measure (IM), which represents the weight that the Transcription Factor has in the prediction of the level of expression of the target gene. We run GENIE3 from the Arboreto library [START_REF] Moerman | GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks[END_REF] using default parameters. GRNBoost2 [START_REF] Moerman | GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks[END_REF] has been developed as a faster alternative to GENIE3. It is thus based on a regression model, using a stochastic gradient boosting machine regression. The inputs and outputs of GRNBoost2 have the same structure of those of GENIE3. Both GRNBoost2 and GENIE3 are part of the SCENIC workflow [START_REF] Aibar | SCENIC: single-cell regulatory network inference and clustering[END_REF]. We run GRNBoost2 from the Arboreto library [START_REF] Moerman | GRNBoost2 and Arboreto: efficient and scalable inference of gene regulatory networks[END_REF] using default parameters.

PPCOR [START_REF] Kim | ppcor: An R Package for a Fast Calculation to Semi-partial Correlation Coefficients[END_REF] infers the presence of a regulatory interaction between two genes by computing the correlation of their expression patterns. To control for possible indirect effects, partial correlation is used instead of a simple correlation, where partial correlation is a measure of the relationship between two variables while controlling for the effect of other variables. The only input of PPCOR is the This is a provisional file, not the final typeset article Running Title expression matrix. The output of PPCOR is a weighted network, where all links are weighted based on the partial correlation between the expression values of the linked nodes/genes. Partial Information Decomposition and Context (PIDC) [START_REF] Chan | Gene Regulatory Network Inference from Single-Cell Data Using Multivariate Information Measures[END_REF] is based on concepts from information theory and uses partial information decomposition (PID) to identify potential regulatory relationships between genes. The only input of PIDC is the expression matrix and its output is a weighted gene-gene network. Context Likelihood of Relatedness (CLR) [START_REF] Faith | Large-Scale Mapping and Validation of Escherichia coli Transcriptional Regulation from a Compendium of Expression Profiles[END_REF] is another commonly-used approach based on concepts from information theory. The measure used by CLR to infer links in between genes is Mutual Information (MI). In contrast with other algorithms also based on MI, such as ARACNE [START_REF] Margolin | ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context[END_REF], CLR adjusts the link weights for the background distribution of the MI values to control for false positives interactions. GeneNet (Opgen-Rhein and Strimmer, 2007) is a method for statistical learning of a high-dimensional causal network. The method first converts a correlation network into a partial correlation graph. Subsequently, a partial ordering of the nodes is established by multiple testing of the log-ratio of standardized partial variances. To make the different network inference algorithms comparable, we applied the same thresholding to all of them, by keeping only the top K links (K=100,000). For methods inferring less than 100,000 links no filtering has been applied.

2.2

Data acquisition and preprocessing Fourteen public scRNA-seq datasets have been used for this benchmark (Table 1): Menon and Lukowski obtained by profiling human retina cells [START_REF] Lukowski | A singlecell transcriptome atlas of the adult human retina[END_REF][START_REF] Menon | Single-cell transcriptomic atlas of the human retina identifies cell types associated with age-related macular degeneration[END_REF]; Zhang and Li profiling T-cells in colorectal cancer (CRC) [START_REF] Li | Reference component analysis of single-cell transcriptomes elucidates cellular heterogeneity in human colorectal tumors[END_REF]Zhang et al., 2019); Hay and Setty profiling human hematopoiesis cells [START_REF] Hay | The Human Cell Atlas bone marrow single-cell interactive web portal[END_REF][START_REF] Setty | Characterization of cell fate probabilities in single-cell data with Palantir[END_REF]. See Table 1 for a complete description of these datasets. The hematopoiesis datasets were split according to their cell type of origin. Only those cell types reported in both studies by Hay et al. And by Setty et al. were considered. We thus obtained a total of ten scRNA-seq datasets in hematopoiesis spanning five cell types: HSC, CLP, Monocyte, Erythroblast and Dendritic Cell. After downloading the data, we filtered the genes based on their total count number (< 3 *0.01*number of cells), as well as on the number of cells in which they are detected (>0.01*number of cells), as described in [START_REF] Aibar | SCENIC: single-cell regulatory network inference and clustering[END_REF]. The gene filtering is performed on each dataset independently. Then, for each biological condition (CRC T-cells, retina, hematopoiesis), only the genes retained for both datasets were selected for network inference. The number of genes retained after filtering are reported in the last column of Table 1. Finally, the data were log2-normalised before applying the different network inference algorithms.

2.3

Indexes employed to measure the reproducibility of the network inference algorithms Percentage of intersection (perINT) and Weighted Jaccard Similarity (WJS) have been employed here to test the reproducibility of the network inference algorithms. The percentage of intersection is used to detect the presence of links shared between two compared networks, while WJS takes into account the similarity of the weights associated with the links shared between the compared networks. Given two networks N1 and N2 inferred respectively from scRNAseq datasets D1 and D2, and indicating as |"| the number of links in the network N, the percentage of intersection (perINT) is computed as: #$%&"'(" ! , " " ) = ) , while the Weighted Jaccard Similarity (WJS) (Tantardini et al., 2019), is defined as

|$ ! ∩$ " | &'((|$ ! |,|$ " |
,-.(" ! , " " ) = ∑ &'(-. # ! ,. $ " / |&| #'! ∑ &01-. # ! ,. $ " / |&| #'!
, where / ! , / " are the vectors of weights associated with the links in common between N1 and N2. In addition, to compare the inferred links to a ground-truth, we considered two additional scores: RcisTarget and Regulatory Circuit scores. We derived the RcisTarget score from the application of the RcisTarget tool [START_REF] Aerts | Robust target gene discovery through transcriptome perturbations and genome-wide enhancer predictions in Drosophila uncovers a regulatory basis for sensory specification[END_REF][START_REF] Aibar | SCENIC: single-cell regulatory network inference and clustering[END_REF]. Given a network of TF-gene interaction, RcisTarget predicts candidate target genes of a TF by looking at the DNA motifs that are significantly over-represented in the surroundings of the Transcription Start Site (TSS) of all the genes that are linked to the TF. We here consider the links validated by RcisTarget as ground-truth and we compare them with the inferred networks, by computing: 0123'4%5$6317%$(" ! ) = (2&345678'(9:;54:4(<∈$ ! ∧?08'@0<4@3ABC':D05E4<

|$ ! | . In the case of the methods inferring links between all genes, a selection of links connecting TFs with possible target genes is performed before computing the RcisTaget score. The Regulatory Circuits score instead is obtained by computing the intersection between an inferred network and tissue-specific regulatory circuits from http://www.regulatorycircuits.org [START_REF] Marbach | Tissuespecific regulatory circuits reveal variable modular perturbations across complex diseases[END_REF]. The regulatory circuits considered are the following: adult retina for retina, lymphocytes for CRC T-cells and CD34 stem cell derived for hematopoiesis. We here computed the Regulatory Circuits score for a network N1 as:

0$589467%:;2%18263317%$(" ! ) = |$ ! ∩'<:0::6C'0<4@54E280<65AC'5C2'<| |$ ! | .

Results

Based on previous works [START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF][START_REF] Pratapa | Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[END_REF], we selected the six single-cell network inference algorithms that do not require an ordering of the cells according to pseudo-time (GENIE3, GRNBoost2, PPCOR, PIDC, CLR and GeneNet see Materials and Methods) and we evaluated them based on their reproducibility, i.e. their ability to infer similar networks once applied to two independent datasets from the same biological condition (e.g. two independent scRNA-seq datasets of colorectal cancer). The reproducibility is measured based on the Percentage of intersection (perINT) and Weighted Jaccard Similarity (WJS) indexes (see Materials and Methods). In addition, we computed the intersection with two instances of ground-truth, based on the RcisTarget and on Regulatory Circuits scores (see Materials and Methods). The evaluation is repeated across three biological conditions: human retina, T-cells in colorectal cancer and human hematopoiesis, for a total of fourteen independent scRNAseq datasets. See Figure 1 for an overview of the benchmark workflow. While in previous benchmarks [START_REF] Chen | Evaluating methods of inferring gene regulatory networks highlights their lack of performance for single cell gene expression data[END_REF][START_REF] Pratapa | Benchmarking algorithms for gene regulatory network inference from single-cell transcriptomic data[END_REF], a low number of highly variable genes had been taken into account (100-1,000 genes), we here tested the ability of the This is a provisional file, not the final typeset article

Running Title

algorithms to infer networks involving all expressed genes (see Materials and Methods for details on the procedure used to filter genes). Indeed, filtering only the top 100-1,000 varying genes is a strong limitation. Restricting the nodes of the inferred network to a low number of genes is reasonable when a manually curated list of relevant genes is available (for example marker genes identified by wet-lab experiments). However, when such a list is not available, working only with the top 100-1,000 varying genes may overlook genes and interactions playing a key role in the regulatory programs of the biological system. We thus tested the various network inference algorithms once applied to scRNAseq datasets containing 6,000-11,000 genes.

In our test cases, PIDC failed to reconstruct networks for two main reasons: (i) the algorithms was slow, especially in the discretization step required to infer a network, and (ii) the use of multivariate information measures impose to have a number of genes much lower than the number of cells, thus requiring to drastically filter out the starting set of genes. Overall, PIDC thus resulted to be more adequate to infer small networks (100-1000 nodes/genes), which are not the focus of this work.

Reproducibility in human retina

We applied GENIE3, GRNBoost2, PPCOR, CLR and GeneNet to two independent scRNA-seq datasets of human retina, reported in [START_REF] Menon | Single-cell transcriptomic atlas of the human retina identifies cell types associated with age-related macular degeneration[END_REF] and in (Lukowsk et al., 2019) (see Materials and Methods). After filtering, the two datasets span 6,212 common genes across a comparable number of cells: 20,091 in Menon versus 20,009 in Lukowski. We thus inferred a total of ten networks. For details on the number of links before and after thresholding refer to Supplementary Table 1. We then evaluated the reproducibility of each algorithm by computing the Percentage of intersection (perINT) and the Weighted Jaccard Similarity (WJS) between the networks inferred independently from the two datasets. While perINT is intended to test the amount of common links between the two networks, the WJS takes also into account the similarity of the weights associated with the common links. As shown in Figure 2A, GENIE3 (45.9% perINT and 0.28 WJS) and GRNBoost2 (41.1% perINT and 0.25 WJS) are the algorithms showing the highest reproducibility, with GENIE3 performing slightly better. At the same time, in agreement with the results of the previous benchmarks, the intersection with the ground truth considered remains rather low but higher for GRNBoost2 (1% RcisTarget score and 4.2% Regulatory Circuits score). Similar performances apply also for the other network inference methods.

Reproducibility in colorectal cancer (CRC) T-cells

We further tested the performances of GENIE3, GRNBoost2, PPCOR, CLR and GeneNet in colorectal cancer (CRC) T-cells. The two datasets used in this case are taken from (Zhang et al., 2019) and [START_REF] Li | Reference component analysis of single-cell transcriptomes elucidates cellular heterogeneity in human colorectal tumors[END_REF] (see Materials and Methods), restricting the last dataset to only T-cells (see Materials and Methods). After filtering, we obtained datasets composed of 11,242 common genes and a widely varying number of cells: 10,805 for Zhang, and 35 for Li. We applied GENIE3, GRNBoost2, PPCOR, CLR and GeneNet independently to the two datasets, for details on the number of links before and after thresholding refer to Supplementary Table 1. Of note, PPCOR has been excluded from this comparison, as it produced partial correlation values outside the range [-1;1] for the Li et al. dataset. After computation of the perINT and WJS (Figure 2B), GENIE3 (3% perINT and 0.008 WJS) and GRNBoost2 (3.4% perINT and 0.007 WJS) emerged as the best performing methods. The reproducibility indexes are quite low in this test case, probably due to the low number of cells present in the Li dataset (35 cells). Also the RcisTarget and Regulatory Circuits scores reflecting the intersection with a ground-truth are quite low for all algorithms, with GRNBoost2 showing better performances (4% RcisTarget score and 14.6% Regulatory Circuits score).

3.3

Reproducibility in human Hematopoiesis Human hematopoiesis has been used as the third biological context for the comparison of GENIE3, GRNBoost2, PPCOR, CLR and GeneNet. The hematopoiesis datasets were split according to the different cell types profiled: HSC, CLP, Monocyte, Erythroblast and Dendritic Cell, obtaining a total of 10 scRNA-seq datasets. Networks were thus inferred on each cell type independently with GENIE3, GRNBoost2, PPCOR, CLR and GeneNet, resulting in a total of 50 networks. Also in this case, details on the number of links before and after thresholding are available in Supplementary Table 1. As for CRC T-cells, PPCOR produced networks composed of links with partial correlation higher than 1 and/or lower than -1 for some CLPs, and Monocytes. For this reason, we did not consider PPCOR in the reproducibility evaluation for these cell types. The reproducibility was then tested for each cell type using the perINT and WJS indexes (Figure 2C-D). GENIE3 displayed the best performances with percentages of intersection of 26-56% and WJS at 0.13-0.37. Consistently with previous observations, the RcisTarget and Regulatory Circuits scores remain low for all cell types and all methods, with GRNBoost2 having slightly better performances than GENIE3 (approx. 2-4.2% and 4-7.6%, respectively) (Figure 2 E,F).

3.4

Stability with respect to link thresholding in the inferred networks In the previous experiments the 100,000 top-ranked links have been taken into account for all methods.

Here we test to which extent our conclusions, regarding the reproducibility of the benchmarked methods, are stable with respect to the number (K) of links retained in each network. We thus apply a more stringent filtering, considered a number (K) of top-ranked links of 10,000, 1000 and 100. GeneNet has been excluded from this analysis as the number of its inferred links is lower than 1000 in most of the cases. After thresholding, the intersection between the networks inferred from independent datasets from the same biological condition were evaluated, using the percINT and WJS as above. As shown in Figure 3, the performances of the various algorithms are quite heterogeneous once different thresholds (K) are considered. As observed in the previous sections, GENIE3 tends to have better performances for high K. However, for low numbers of links (K=1000-100) GRNBoost2 and CLR tend to predominate in most of the cases.

Stability with respect to technical variations in the input data: number of profiled cells, sequencing platform and cell type annotation

In the experiments performed above, we tested the reproducibility of the network inference algorithms by using two independent datasets for each biological condition (e.g. human retina). A limitation of this approach comes from the technical differences between the protocols followed to generate these datasets: different sequencing platforms, different procedures used for the annotation of the cell types, and different number of cells. All these technical differences could impact our results. To evaluate the stability of the results against technical variations, we used the largest dataset, from [START_REF] Menon | Single-cell transcriptomic atlas of the human retina identifies cell types associated with age-related macular degeneration[END_REF], encompassing 20,091 cells. We splitted this dataset into two subsets, keeping the proportions of the various cell types constant. We then applied the five network inference algorithms independently to the two subsets, and we evaluated the reproducibility of the algorithms using perINT and WJS, as in the previous tests. To further assess the effect of the number of cells on network This is a provisional file, not the final typeset article Running Title inference, we split the same scRNAseq dataset generated by Menon et al. three times to obtain couples of datasets encompassing decreasing number of cells: 10,000, 1,000 and 100. Note that for all these comparisons, the sequencing platform and/or the method/technique used to annotate the cells are identical for all subsets. PPCOR inferred networks for 10,000 and 1,000 cells, but failed at 100 cells (see Supplementary Table 2). For details on the number of links before and after thresholding (K=100,000) refer to Supplementary Table 2. Overall, as shown in Figure 4, GENIE 3 emerged again as the best performing method in all cases. Of note, for low number of cells, a general decrease in reproducibility is observed for all network inference methods, which can be justified by a lower accuracy in the link estimation due to the low number of observations (cells).

The scNET Jupyter notebook

To foster the reproducibility of all the results and figures presented in this study, we implemented the corresponding code in a Jupyter notebook, which is available at https://github.com/ComputationalSystemsBiology/scNET, together with a Conda package containing all the required libraries. Importantly, scNET can be used to benchmark new network inference algorithms based on their reproducibility, or further test GENIE3, PPCOR, GRNBoost2, CLR and GeneNet on user-provided datasets. 4 Discussion Starting from the benchmark of Paratapa et al. ( 2020), we evaluated the network inference algorithms from a complementary perspective by assessing their reproducibility. We were thus interested to test if the algorithms would infer similar networks when applied to pairs of independent datasets from the same biological condition (e.g. T-cells in colorectal cancer). Our benchmark focused on real patientderived data spanning three biological contexts: human retina, T-cells in CRC, and human hematopoiesis cells. We thus consider highly different biological contexts, going from cancer tissue, to isolated healthy immune cells, and to a mixture of normal retina cells combined in a single dataset. Importantly, we aimed at inferring networks involving a much higher number of genes compared to previous works. In agreement with previous benchmarks, all network inference algorithms generated networks having low intersections with ground-truth. Of note the ground-truth considered here, based on RcisTarget and regulatory circuits, is different and complementary to those used in previous benchmarks. This disappointing result might arise for different reasons, potentially adding up. Limitations can be present in the input data, as scRNAseq may not provide sufficient resolution for reliable network inference, and technical and experimental factors present in the input data might affect information content. Turning to the inference algorithm, limitations may arise from underlying statistical assumptions and the documented lack of uniqueness in the solution of the network inference problem. Finally, the ground-truth network considered here and in previous benchmarks may not be sufficiently comprehensive. PPCOR provided weights outside the normal range of correlation values ([-1,1]) for datasets having less than 1,000 cells. Such inconsistencies are likely due to numerical problems arising when the input dataset encompasses many more genes than cells. PIDC was the algorithm that suffered the most when applied to high numbers of genes. Overall, for high link numbers (K=100,000) GENIE3 consistently generated the most reproducible results across all the three biological contexts considered. Furthermore, its performances proved to be stable with respect to the single-cell sequencing platform, the cell type annotation system, the number of cells considered. Once instead a more stringent filtering is considered (K=100-100) CLR and GRNBoost2 show better performances. However, even the best performing methods show reproducibility scores that are less than ideal (26-54% perINT and 0.1-0.3 WJS), indicating that lot of improvement is still needed in the design of network inference methods for scRNA-seq data. We here considered network inference methods that are highly heterogeneous. With some algorithms, as PPCOR and GeneNet, inferring links between all possible couples of genes and others, as GENIE3 and GRNBoost2, only inferring links between Transcription Factors (TFs) and possible target genes. We tried to make the inferred networks comparable by fixing the number of links in all network to a certain value K, thus obtaining networks with the same density. However, in principle, methods inferring only TF-target links will have higher chances to be reproducible in our comparison. At the same time, once the links of PPCOR and GeneNet are restricted to only TF-target links, the dimension of the networks drastically decreases (sometimes empty networks are obtained) and their performances do not improve. The main limitation of this benchmark is the number of considered network inference algorithms. Future extensions of this study could include pseudotime-based network inference methods, once adequate datasets will become available. To date, available independent datasets relevant for pseudotime-based network inference algorithms (e.g. cells profiled during development stimulation) present too many experimental variations to be employed for a reliable evaluation of reproducibility. Of note, such extensions will be greatly facilitated by taking advantage of the Jupyter notebook (scNET) provided as supplementary material. This is a provisional file, not the final typeset article This is a provisional file, not the final typeset article

Running Title

Sonawane, A. R., Weiss, S. T., Glass, K., and Sharma, A. (2019)

Figure 2 .

 2 Figure 2. Reproducibility performances of the various network inference algorithms across the three biological contexts: human retina, colorectal cancer T-cells and human hematopoiesis. Panels A and B report the Percentage of intersection (perINT), Weighted Jaccard Similarity (WJS), RcisTarget score and Regulatory Circuits score obtained by each the benchmarked algorithm (GRNBoost2, GENIE3, PPCOR, CLR and GeneNet) in human retina and colorectal cancer T-cells, respectively. Panels C-F summarize the performances of the same algorithms in hematopoiesis, with perINT (C), WJS (D), RcisTarget score (E) and Regulatory Circuits score (F).Figure 3. perINT and WJS according to different network thresholding. The perINT and WJS are reported for varying thresholds on the weight distribution of the links of the inferred networks. The results are reported for all the tested datasets (A) retina, (B) CRC T-cells, (C ) CLPs, (D) Dendritic cells, (E) Erythrocytes, (F) HSCs, (G) Monocytes.Figure 4. Stability of the network inference performances with respect to technical variations in the input data. Reproducibility scores of GRENBoost2 (red), GENIE3 (black) PPCOR (yellow), CLR (gray) and GeneNet (green) across different splittings of the Menon retina dataset. Panels A and B report the percentage of intersection (perINT) and Weighted Jaccard Similarity (WJS), respectively.

Figure 3 .

 3 Figure 2. Reproducibility performances of the various network inference algorithms across the three biological contexts: human retina, colorectal cancer T-cells and human hematopoiesis. Panels A and B report the Percentage of intersection (perINT), Weighted Jaccard Similarity (WJS), RcisTarget score and Regulatory Circuits score obtained by each the benchmarked algorithm (GRNBoost2, GENIE3, PPCOR, CLR and GeneNet) in human retina and colorectal cancer T-cells, respectively. Panels C-F summarize the performances of the same algorithms in hematopoiesis, with perINT (C), WJS (D), RcisTarget score (E) and Regulatory Circuits score (F).Figure 3. perINT and WJS according to different network thresholding. The perINT and WJS are reported for varying thresholds on the weight distribution of the links of the inferred networks. The results are reported for all the tested datasets (A) retina, (B) CRC T-cells, (C ) CLPs, (D) Dendritic cells, (E) Erythrocytes, (F) HSCs, (G) Monocytes.Figure 4. Stability of the network inference performances with respect to technical variations in the input data. Reproducibility scores of GRENBoost2 (red), GENIE3 (black) PPCOR (yellow), CLR (gray) and GeneNet (green) across different splittings of the Menon retina dataset. Panels A and B report the percentage of intersection (perINT) and Weighted Jaccard Similarity (WJS), respectively.

Figure 4 .

 4 Figure 2. Reproducibility performances of the various network inference algorithms across the three biological contexts: human retina, colorectal cancer T-cells and human hematopoiesis. Panels A and B report the Percentage of intersection (perINT), Weighted Jaccard Similarity (WJS), RcisTarget score and Regulatory Circuits score obtained by each the benchmarked algorithm (GRNBoost2, GENIE3, PPCOR, CLR and GeneNet) in human retina and colorectal cancer T-cells, respectively. Panels C-F summarize the performances of the same algorithms in hematopoiesis, with perINT (C), WJS (D), RcisTarget score (E) and Regulatory Circuits score (F).Figure 3. perINT and WJS according to different network thresholding. The perINT and WJS are reported for varying thresholds on the weight distribution of the links of the inferred networks. The results are reported for all the tested datasets (A) retina, (B) CRC T-cells, (C ) CLPs, (D) Dendritic cells, (E) Erythrocytes, (F) HSCs, (G) Monocytes.Figure 4. Stability of the network inference performances with respect to technical variations in the input data. Reproducibility scores of GRENBoost2 (red), GENIE3 (black) PPCOR (yellow), CLR (gray) and GeneNet (green) across different splittings of the Menon retina dataset. Panels A and B report the percentage of intersection (perINT) and Weighted Jaccard Similarity (WJS), respectively.

Datasets employed in this benchmark Data Name Biological context Sequencing technology Number of cells Cell type annotation strategy Associated publication Number of genes after preprocessing

  . Network Medicine in the Age of Biomedical Big Data. Frontiers in Genetics 10, 294. doi:10.3389/fgene.2019.00294. Tantardini, M., Ieva, F., Tajoli, L., and Piccardi, C. (2019). Comparing methods for comparing networks. Scientific Reports 9, 17557. doi:10.1038/s41598-019-53708-y.

	The DREAM5 Consortium, Marbach, D., Costello, J. C., Küffner, R., Vega, N. M., Prill, R. J., et al.
		(2012). Wisdom of crowds for robust gene network inference. Nat Methods 9, 796-804.
		doi:10.1038/nmeth.2016.
	Verny, L., Sella, N., Affeldt, S., Singh, P. P., and Isambert, H. (2017). Learning causal networks with
		latent variables from multivariate information in genomic data. PLoS Comput Biol 13,
		e1005662. doi:10.1371/journal.pcbi.1005662.
	Zhang, Y., Zheng, L., Zhang, L., Hu, X., Ren, X., and Zhang, Z. (2019). Deep single-cell RNA
		sequencing data of individual T cells from treatment-naïve colorectal cancer patients. Sci
		Data 6, 131. doi:10.1038/s41597-019-0131-5.
	Tables
	Table1. Menon Human retina	10X Genomics	20,091	manually curated	(Menon et al.,	6212
					marker genes	2019)
	Lukowski	Human retina	10X Genomics	20,009	no annotation	(Lukowski et	6212
						al., 2019)
	Zhang	CRC T-cells	Smart-Seq2,	10,805	FACS sorted	(Zhang et al.,	11242
						2019)
	Li	CRC T-cells	HiSeq 2000	375 cells (of	manually curated	(Li et al., 2017)	11242
			Illumina	which 35 T-	marker genes
				cells)
	Hay	human	10X Genomics	101,935	MarkerFinder ICGS	(Hay et al.,	7038
	hematopoiesis				2018)
	Setty	human	10X Genomics	12,046	Sorted bulk	(Setty et al.,	7038
	hematopoiesis			hematopoietic	2019)
					populations
	Figures Legends
	Figure 1.

Summary of the workflow followed in this benchmark.

Acknowledgments

We thank the bioinformatics platform of IBENS for the computational/infrastructural support. We thank Michael Mason, Anaïs Baudot and Sabine Tejpar for the scientific feedbacks on the work.

Data Availability Statement

The datasets for this study can be accessed from their associated publications (see Table1). All the analyses are reproducible using the scNET Jupyter notebook available at https://github.com/ComputationalSystemsBiology/scNET.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Author Contributions

LC designed the analysis. YK performed the analysis. LC and DT co-supervised the study. All authors contributed to the manuscript and approved the submitted version.


MD Claire Hentzen

Anaïs Villaumé

MD Nicolas Turmel

MD Camille Chesnel

MD Frédérique Le Breton

MD Rebecca Haddad

PhD Gérard Amarenco

Are falls in people with multiple sclerosis related to the severity of urinary disorders?

Keywords: fall, multiple sclerosis, lower urinary tract symptoms, overactive bladder

Background. Lower urinary tract symptoms (LUTSs), especially overactive bladder, are frequent in people with multiple sclerosis (PwMS). Urinary urgency and urge urinary incontinence could lead to precipitation and thus could increase the risk of falling in these individuals.

Objective. We aimed to assess the association between severity of LUTSs and risk of falling in PwMS.

Methods

. PwMS with LUTSs were recruited in a neuro-urology department. Participants were asked about the number of falls in the past 3 months and their circumstances. Severity of LUTSs was assessed by the Urinary Symptoms Profile (USP) score, and individuals were classified as with or without urinary incontinence.

). Overall, 20 (13%) patients reported one fall during the past 3 months, and 43 (28%) reported at least 2 falls. Only 9 (6%) patients reported a fall on the way to the toilet, 6 during a urinary urgency. No link was found between falls and urinary incontinence (p=0.71), type or severity of urinary symptoms (overactivity, voiding dysfunction or stress incontinence, p > 0.05). Falls on the way to the bathroom was associated with high USP score related to overactive bladder (p=0.03) and severe nocturia (> 2 nocturnal micturitions) (p<0.01). Falling at night was also associated with severe nocturia (p < 0.001).

Conclusions. The severity of LUTSs and presence of urinary incontinence do not appear related to the risk of falling in PwMS and urinary disorders but rather to the specific risk of falling on the way to the bathroom. Severe nocturia increases the risk of falling at night. Further studies are needed to assess the impact of LUTS treatment on the risk of falling. ClinicalTrials.gov (NCT04338646)

Introduction

Falls are a common problem in people with multiple sclerosis (PwMS) and can lead to severe consequences (trauma, fear of falling, reduction of social activities). Prevention of falls is one of the priority targets of rehabilitation for PwMS and walking difficulties, which can result from different factors (motor impairment, ataxia, visual alterations, sensitive disorders, fatigability etc.). Some clinical characteristics related to motor impairment (gait velocity, 12-Item Multiple Sclerosis Walking Scale [MSWS12] score) [START_REF] Fritz | Distinguishing among multiple sclerosis fallers, near-fallers and non-fallers[END_REF] or fear of falling (Falls Efficacy Scale-International [FES-I] score) [START_REF] Gunn | Relationship Between Physiological and Perceived Fall Risk in People With Multiple Sclerosis: Implications for Assessment and Management[END_REF] have been identified as risk factors of falling. Urinary incontinence has also been considered to predict falls, initially in older adults [START_REF] Chen | Development of a fall-risk assessment profile for community-dwelling older adults by using the National Health Interview Survey in Taiwan[END_REF][START_REF] De Rekeneire | Is a fall just a fall: correlates of falling in healthy older persons. The Health, Aging and Body Composition Study[END_REF] but secondary in PwMS [START_REF] Matsuda | Falls in multiple sclerosis[END_REF][START_REF] Sung | Bladder function and falls in individuals with multiple sclerosis[END_REF][START_REF] Gunn | Risk factors for falls in multiple sclerosis: an observational study[END_REF].

However, lower urinary tract symptoms (LUTSs) are frequent in PwMS, with a high prevalence (32-96.8%) and an increase with MS duration and severity of neurological deficiencies and disabilities [START_REF] De Sèze | The neurogenic bladder in multiple sclerosis: review of the literature and proposal of management guidelines[END_REF]. Overactive bladder with urgency, frequency, urge urinary incontinence and nocturia is the most common symptom, reported by 37% to 99% of PwMS. Overactive bladder affects quality of life, and a neurogenic bladder with detrusor overactivity associated with detrusor sphincter dyssynergia on urodynamic assessment can lead to an altered upper urinary tract (reflux, dilatation, urinary tract infection) [START_REF] Giannantoni | Urological dysfunctions and upper urinary tract involvement in multiple sclerosis patients[END_REF][START_REF] Onal | Voiding dysfunction due to multiple sclerosis: a large scale retrospective analysis[END_REF][START_REF] Krhut | Analysis of the upper urinary tract function in multiple sclerosis patients[END_REF][START_REF] Castel-Lacanal | Urinary complications and risk factors in symptomatic multiple sclerosis patients. Study of a cohort of 328 patients[END_REF][START_REF] Shakir | Assessment of Renal Deterioration and Associated Risk Factors in Patients With Multiple Sclerosis[END_REF]. Despite its high prevalence and impact on quality of life, the severity of LUTSs has never been studied as specific risk factor of falling.

However, urinary urgency and urinary incontinence could lead to precipitation and thus could increase the risk of falling in these patients. This hypothesis had been previously evoked in older adults [START_REF] Paquin | The effects of a strong desire to void on gait for incontinent and continent older communitydwelling women at risk of falls[END_REF][START_REF] Roggeman | The role of lower urinary tract symptoms in fall risk assessment tools in hospitals: a review[END_REF].

The primary aim of the study was to assess the association between severity of LUTSs and risk of falling in PwMS. The secondary aim was to correlate patient characteristics and clinical tests with risk of falling. The study was conducted and results were reported according to the STROBE guidelines for cross-sectional studies (see additional material).

Material and methods

Study design

This study was approved by the local ethics review board (RCB: 2018-A01644-51) and all participants provided written informed consent before inclusion in this cross-sectional study.

This study was registered at ClinicalTrials.gov (NCT04338646). Participants were recruited in a neuro-urology department during a medical appointment or a urodynamic assessment or in a day hospital related to urinary disorders between July 2019 and January 2020.

Participants

Inclusion criteria were age ≥ 18 years with an MS diagnosis, LUTSs with or without treatment, and Expanded Disability Status Scale (EDSS) score 1 to 6.5. Exclusion criteria were relapse of MS in the past month, urinary tract infection on the day of inclusion, and inability to complete the symptom questionnaires.

Procedure

Collected data were age, sex, EDSS score, body mass index, need for visual correction and use of a walking aid. Patients were asked about the number of falls in the past 3 months and in the past year and the circumstances in which they occurred (frequency, home, outdoors, going to void, during urinary urgency, nocturia). Number of falls was classified as 0, 1 or at least 2 to quantify the severity of risk of falling. A better precision was not possible owing to the retrospective collection and would have led to classification bias. Severity of LUTSs was assessed by the Urinary Symptoms Profile (USP) score [START_REF] Haab | Comprehensive evaluation of bladder and urethral dysfunction symptoms: development and psychometric validation of the Urinary Symptom Profile (USP) questionnaire[END_REF] and participants were classified as with or without urinary incontinence. Number of micturitions at night was specifically asked.

To take into account motor difficulties and fear of falling, other clinical evaluations were done.

The impact of MS on walking was assessed with the MSWS12 questionnaire [START_REF] Hobart | Measuring the impact of MS on walking ability: the 12-Item MS Walking Scale (MSWS-12)[END_REF], the EDSS score, and a clinical test with the Time to be Ready to Void (TRV) test [START_REF] Hentzen | Time to be Ready to Void: A new tool to assess the time needed to perform micturition for patients with multiple sclerosis[END_REF]. Fear of falling was assessed by a simple question and with the FES-I questionnaire [START_REF] Van Vliet | Falls efficacy scaleinternational: a cross-sectional validation in people with multiple sclerosis[END_REF].

The primary aim was to assess the association between severity of LUTSs and occurrence of falls during the past 3 months. The primary outcome was the importance of overactive bladder (OAB) symptoms with the OAB USP score. Secondary outcomes were the existence of urinary incontinence, the warning time (defined as the time from the first sensation of urgency to voiding or incontinence), the importance of nocturia and the other sub-scores of the USP questionnaire (voiding dysfunction and stress urinary incontinence).

Secondary aims were to look for the relationship between severity of LUTSs and specific circumstances of falls, severity of LUTS and occurrence of falls during the past year, and to assess the association between falls and classical risk factors of falls.

Statistical analysis

Because the frequency of falls in PwMS is variable in literature (30-63% [START_REF] Giannì | A systematic review of factors associated with accidental falls in people with multiple sclerosis: a meta-analytic approach[END_REF]), and there are no data on the impact of the severity of urinary disorders on the risk of falling, the number of participants to be included could not be calculated. The minimal target of inclusion was 100 participants by estimating the reported fall rate at about 30%, which allowed us to use parametric tests on the primary outcome. The number of missing data was very low, so these were not replaced. Statistical analyses involved using R for Windows (Rx64

Results

Patient characteristics

In total, 509 PwMS were seen in the department during the inclusion period; 167 did not meet the inclusion criteria and 188 declined participation. We included 154 patients (69% women; Figure), with mean (SD) age 50.1 (11.5) years and median EDSS 5 (range 3-6). All participants completed the full study. The demographic characteristics of participants are reported in Table 1. Over the past 3 months, 13% (n=20) of PwMS reported 1 fall and 28% (n=43) reported at least 2 falls; 25% reported falls indoors and 30% outdoors (Table 2). Over the past year, 10% reported 1 fall and 45% reported at least 2 falls; 36% reported falls indoors and 43% outdoors.

Association between reported falls and LUTSs

We found no link between the reported falls during the past 3 months and severity of LUTSs (USP sub-scores: OAB, voiding dysfunction or stress incontinence, p > 0.05) or urinary incontinence (p=0.71). Falling at night was associated with severe nocturia (> 2 nocturnal micturitions) (p < 0.001) but not with overall risk of falling (including falls during the day and at night). Similar results were found between the reported falls during the past year and severity of LUTSs. Results are summarized in Table 3.

Association between reported falls on the way to the bathroom and LUTSs

Over the past 3 months, 9 (6%) patients reported a fall on the way to the bathroom, 6 of them during a urinary urgency. These 9 patients had higher mean (SD) OAB USP score than others (9.9 [5.8] vs 6.5 [4.7], p = 0.03) and 56% (n=5) had severe nocturia (vs 16% for the rest of the participants, p <0.01). We found no difference concerning urinary incontinence or use of selfcatheterization.

Over the past year, 11 (7%) patients reported at least one fall on the way to the bathroom, 7 of them during a urinary urgency, which represent 13% of participants who had ever fallen. Again, participants who reported fall on the way to the bathroom had higher mean OAB USP score than others (11.3 [3.8] vs 6.5 [4.7], p < 0.01) and a higher proportion had severe nocturia (45% vs 17%, p = 0.04).

Association between reported falls and motor symptoms or fear of falling

On univariate analysis, factors associated with increased risk of falling during the past 3 months were EDSS score (p<0.001), fear of falling and FES-I score (p<0.001), MSWS12 score (p<0.0001), use of a walking device (p=0.02), gait velocity (mobility stage of TRV) (p<0.01), settled stage of TRV (p<0.0001), and total TRV score (p<0.001). The same factors were significantly associated with the risk of fall during the past year.

On multivariate analysis, falls in the past 3 months was associated with MSWS12 score (odds ratio [OR] 1.04 [95% confidence interval (CI) 1.01-1.06], p < 0.001) and total TRV score (1.03 [1.0-1.06], p = 0.03) (Table 4). Falls in the past year was associated with fear of falling (3.15

Discussion

Among the 154 PwMS with urinary disorders, 63 (41%) reported at least one fall in the past 3 months, and 9 of these 63 (14%) patients reported at least one fall while they were going to the bathroom. Severity of LUTSs did not seem to be linked with risk of falls in general but was linked to the specific risk of fall on the way to the bathroom. Severe nocturia was associated with risk of falls during the night. High MSWS12 and TRV scores were associated with falls in the past 3 months on multivariate analysis.

Previous studies have explored the frequency of falls in PwMS, with or without urinary disorders. A systematic review reported a proportion of fallers between 30% and 63% in a pool of 2,425 patients [START_REF] Giannì | A systematic review of factors associated with accidental falls in people with multiple sclerosis: a meta-analytic approach[END_REF]. Most falls were reported retrospectively by the patients, over a period from 1 to 12 months depending on the study. Only 3 of 15 studies were prospective. In comparison with large prospective cohorts, we found a slightly lower rate of falls. Thus, in a cohort of 416 PwMS with median EDSS close to 4, Gunn et al. reported 17% of patients with one fall and 37% with at least 2 falls in 3 months [START_REF] Gunn | Relationship Between Physiological and Perceived Fall Risk in People With Multiple Sclerosis: Implications for Assessment and Management[END_REF]. In another cohort of 537 PwMS, with a median EDSS similar to the one in our study, 56% patients reported falls, and 37% mentioned at least 2 falls during a 3-month period [START_REF] Nilsagård | Falls in people with MS--an individual data meta-analysis from studies from Australia, Sweden, United Kingdom and the United States[END_REF]. Most falls were indoors. Regarding falls over the past year, Kalron et al. found 47% of patients reporting at least 2 falls, close to our rate of 45% [START_REF] Kalron | The relationship between specific cognitive domains, fear of falling, and falls in people with multiple sclerosis[END_REF]. However, the mean EDSS was lower in their study. Several hypotheses may explain these differences: a memory bias in the number of falls reported retrospectively (underestimation) or possible differences in MS typology (duration, course, clinical presentation etc.).

Concerning the impact of LUTS on the risk of falling, our results are not completely consistent with the current literature. In a study of 92 individuals with MS, Sung et al. reported an increased fall frequency for patients with severe urinary disorders versus those with mild troubles on multivariate analysis [START_REF] Sung | Bladder function and falls in individuals with multiple sclerosis[END_REF]. However, the severity of urinary disorders was assessed with the EDSS sub-item related urinary continence, which is not precise or specific and could explain our contradictory results. Other studies found a link between risk of falling and urinary incontinence, but not all patients had LUTSs [START_REF] Finlayson | Risk factors for falling among people aged 45 to 90 years with multiple sclerosis[END_REF], and the results were not adjusted by MS severity or disability [START_REF] Fritz | Distinguishing among multiple sclerosis fallers, near-fallers and non-fallers[END_REF][START_REF] Matsuda | Falls in multiple sclerosis[END_REF]. In older adults, urinary incontinence is an independent risk factor for falls [START_REF] De Rekeneire | Is a fall just a fall: correlates of falling in healthy older persons. The Health, Aging and Body Composition Study[END_REF]. A systematic review reported an odds of falling of 1.54 (95% CI 1.41-1.69) for urge incontinence and 1.11 (1.00-1.23) for stress incontinence [START_REF] Chiarelli | Urinary incontinence is associated with an increase in falls: a systematic review[END_REF]. Age is classically also associated with falls [START_REF] Chen | Development of a fall-risk assessment profile for community-dwelling older adults by using the National Health Interview Survey in Taiwan[END_REF]. We did not find an impact of age, probably because our population was young, and severity of MS is the most important factor related to the risk of fall. Finally, some patients had already received treatment for LUTSs, which could reduce urinary symptoms, principally urgency and urge incontinence, and thus, the necessity of precipitation to go to the bathroom.

The severity of LUTSs depends on treatment efficacy and this can induce a bias in symptoms' evaluation. However, treatment of LUTSs was not taken into account in the previous studies either.

In our study, nocturia was associated with falls on the way to the bathroom. The impact on falls under any circumstances may have been underestimated owing to lack of power or to the predominant role of gait disorders in our population. One possible hypothesis is that nocturia may affect sleep and induce increased fatigue, which may contribute to the overall risk of falling throughout the day. Our study found an association between nocturia and risk of falling at night but not the risk of falling in general, but in older adults, nocturia is a well-established risk factor of falls, also associated with an over-risk of fracture [START_REF] Pesonen | The Impact of Nocturia on Falls and Fractures: A Systematic Review and Meta-Analysis[END_REF].

Our study is the first to specifically assess whether the severity of LUTSs increases the risk of falls in patients with MS and urinary disorders. Our population was more homogenous than in previous studies because all participants had LUTSs, which reinforces our results. Urinary symptoms were for the first time assessed with a validated questionnaire, in addition to a simple declarative response of the presence of urinary incontinence. This is also the first time that the circumstances of falls specifically linked with urinary disorders were studied in this population.

Only 9 PwMS reported falls when going to the bathroom and 6 when having a urinary urgency, representing 23% and 15%, respectively, of PwMS who reported falls indoor. In older adults, a fall on the way to the bathroom is a classical circumstance, especially during hospitalization, estimated at 20% to 50% [START_REF] Roggeman | The role of lower urinary tract symptoms in fall risk assessment tools in hospitals: a review[END_REF]. In Parkinson's disease, Sakushima et al. showed that 13.8% of falls occurred on the way to the toilet [START_REF] Sakushima | Influence of urinary urgency and other urinary disturbances on falls in Parkinson's disease[END_REF]. They also found an association with urinary urgency and falls. Thus, going to the bathroom is probably a prevalent circumstance of falls indoor, and the risk could be increased in case of frequency and urge urinary incontinence.

We found high MSWS12 and TRV scores associated with falls in the past 3 months. The MSWS12 score is a well-known prognostic factor for the risk of falling [START_REF] Kalron | The relationship between specific cognitive domains, fear of falling, and falls in people with multiple sclerosis[END_REF][START_REF] Tajali | Predicting falls among patients with multiple sclerosis: Comparison of patient-reported outcomes and performance-based measures of lower extremity functions[END_REF]. This selfreporting questionnaire reports patient perception of gait difficulties. In our study, MSWS12 score was associated with reported falls with an OR close to 1. However, the association was found for both falls in the past 3 months and in the past year, which strengthens our results.

The TRV test is a new validated test for PwMS [START_REF] Hentzen | Time to be Ready to Void: A new tool to assess the time needed to perform micturition for patients with multiple sclerosis[END_REF]. This test assesses quantitatively the time required to reach the bathroom and be ready to void, a complex task, including stand up for a chair, walking, opening a door, and undressing. Several functions are implicated: gait, dexterity, executive functions etc. Because this task is more complex than a Timed up and Go test and because a study showed that the Timed up and Go test failed to predict falls in PwMS [START_REF] Quinn | Discriminative ability and clinical utility of the Timed Up and Go (TUG) in identifying falls risk in people with multiple sclerosis: a prospective cohort study[END_REF], we decided to use the TRV. The ORs for the 2 factors linked with falls were low, probably because of interactions between them and other factors such as EDSS and fear of falling.

Several limitations exist. First, a number of patients declined participation for different reasons (personal or professional imperatives, tiredness secondary to the appointment, did not want to participate), and we cannot exclude a selection bias. Also, the number and circumstances of falls were declarative. This can induce a memory bias, especially in patients with MS and possible cognitive disorders. Moreover, the number of falls during the night or during urinary urgency being low must be taken into account in the interpretation of our results. Second, the influence of urinary treatment was not evaluated. However, the treatment of LUTSs was not modified in the past 3 months, which limits its impact. Therefore, the number of falls in the past 3 months were primarily studied, limiting the influence of treatment, variability of urinary symptoms and modification in motor impairment (e.g., owing to relapse of MS). Finally, because of the heterogeneity of PwMS in neurological presentation, the detailed data from neurological examinations were not recorded. Thus, we could not assess the risk of falling function of some clinical forms (e.g., pyramidal, cerebellar) and the results of this study must be considered by taking into account these different biases. Other factors may come into play, such as the attentional impact of the need to urinate, which in other populations such as older people or those with Parkinson's disease could be important, so these results cannot be generalized to other populations at risk of falling. Many perspectives exist. It could be interesting to assess the impact of the treatment of LUTSs on risk of falling with a prospective study. The impact of increased disability should be taken into account with progressive MS or new relapse. The correlation with manometric results, such as the existence of detrusor overactivity on cystometry, could also be assessed. A specific study of individuals with high EDSS score could be interesting, because these people are more at risk of falling, and some treatments can improve LUTSs. To demonstrate a benefit of treatment of OAB symptoms may lead neurologists to detect these symptoms and address the patient to the neuro-urologist. Finally, a prospective study using the TRV test could be interesting to assess its ability to predict falls in MS.

Conclusion

Severity of LUTSs and urinary incontinence do not seem related to risk of falling in general in PwMS with urinary disorders but are related to increased risk of falling in specific circumstances. Thus, falling on the way to the bathroom is associated with severe OAB symptoms, and severe nocturia increases the risk of falling at night. High MSWS12 and TRV scores were related to falls. Further studies are needed to assess the impact of LUTS treatment on the risk of falling in PwMS.

Legends

[ 1 .

 1 15-8.31], p=0.03) and MSWS12 score (1.03 [1.01-1.05], p < 0.01).

Figure .

 . Figure. Flow chart of the study. The exact reasons for refusing to participate were not

  The impact of the different factors on the risk of falling was assessed by Student t test and ANOVA for quantitative variables if normally distributed or Welch's t-test otherwise, and chi-square test or Fisher exact test for categorical variables. Statistical significance was set at p < 0.05. A logistic regression model was applied to search for risk factors in multivariate analysis.

	3.4.2, R
	Foundation for Statistical Computing, Vienna, Austria). Descriptive data are presented as
	number (%) or mean (SD) for continuous data and median (range) for ordinal data and data not
	normally distributed.

Table 1 .

 1 Initial characteristics of participants with multiple sclerosis.

	Age (years), mean (SD)	50.1 (11.5)
	Female	107 (69%)
	BMI, mean (SD)	24.6 (5.9)
	EDSS, median [Q1-Q3]	5 [3-6]
	MSWS12, mean (SD)	65.2 (22.5)
	FES-I, mean (SD)	31.7 (10.7)
	TRV total (sec), mean (SD)	29.7 (18.8)
	TRV mobility stage (sec), mean (SD)	11.1 (10.2)
	USP score, mean (SD)
	USP stress score (/9)	1 (2)
	USP OAB score (/21)	6.7 (4.7)
	USP voiding dysfunction score (/9)	5 (3.4)
	Warning time (min)
	 1	36 (23%)
	1-5	49 (32%)
	6-15	23 (15%)
	> 15	46 (30%)
	Micturition status
	Spontaneous void	94 (61%)
	CISC	31 (20%)
	Mixed	29 (19%)
	Urinary incontinence
	Indoor	72 (47%)
	Outdoor	88 (57%)
	Nocturia
	0-1	89 (58%)
	2	34 (22%)
	> 2	31 (20%)
	BMI, body mass index; EDSS, Expanded Disability Status Scale; MSWS12, 12-Item Multiple
	Sclerosis Walking Scale; FES-I, Falls Efficacy Scale-International questionnaire; TRV, Time
	to be Ready to Void; USP, Urinary Symptom Profile; OAB, overactive bladder; CISC, Clean
	Intermittent Self-Catheterization

Table 2 .

 2 Falls and their circumstances in the past 3 months and the past year

	Past 3 months	Past year

Table 3 .

 3 Association between falls and lower urinary tract symptoms. Data are mean (SD) unless indicated. USP, Urinary Symptom Profile; OAB, Overactive Bladder; SUI, Stress Urinary Incontinence a ANOVA; b Fisher exact test

		3 months				1 year
		0	1	2	p	0	1	2	p
	USP OAB	6.89 (4.62) 6.0 (5.2)	6.67 (4.66) 0.75 a 6.43 (4.58) 6.47 (5.58) 7.04 (4.64) 0.73 a
	USP SUI	0.96 (2.01) 1.35 (2.3)	0.95 (2.07) 0.73 a 0.84 (2.0)	0.93 (1.98) 1.19 (2.14) 0.61 a
	USP voiding	4.59 (3.45) 4.65 (3.45) 5.88 (3.19) 0.12 a 4.42 (3.44) 4.47 (3.58) 5.60 (3.28) 0.11 a
	dysfunction
	Urinary incontinence
	Yes	45 (49%)	9 (45%)	18 (42%)	0.71 b 31 (45%)	6 (40%)	35 (50%)	0.72 b
	No	46 (51%)	11 (55%)	25 (58%)		38 (55%)	9 (60%)	35 (50%)
	Warning time (min)
	< 1	23 (25%)	4 (20%)	9 (21%)	0.75 b 15 (22%)	4 (27%)	17 (24%)	0.60 b
	1-5	30 (33%)	5 (25%)	14 (32%)		25 (36%)	3 (20%)	21 (30%)
	6-15	10 (11%)	4 (20%)	9 (21%)		7 (10%)	2 (13%)	14 (20%)
	> 15	28 (31%)	7 (35%)	11 (26%)		22 (32%)	6 (40%)	18 (26%)
	Nocturia  2
	Yes	16 (18%)	4 (20%)	10 (23%)	0.74 b 12 (17%)	4 (27%)	14 (20%)	0.71 b
	No	75 (82%)	16 (80%)	33 (77%)		57 (83%)	11 (73%)	56 (80%)

Table 4 .

 4 Association between falls and clinical factors. Data are mean (SD) unless indicated. EDSS, Expanded Disability Status Scale; FES-I, Falls Efficacy Scale-International; MSWS12, 12-Item Multiple Sclerosis Walking Scale; TRV, Time to be Ready to Void a Student t-test b Welch two sample t-test; c, chi-squared test; LR, logistic regression (multivariate analysis); OR, odds ratio; 95% CI, 95% confidence interval

	3 months


Guillaume Forbin
email: g.forbin@surrey.ac.uk-bbessere@univ-lr.fr-boldys@utia.cas.cz-david.tschumperle@greyc.ensicaen.fr

Bernard Besserer

Jiří Boldyš

David Tschumperlé

TEMPORAL EXTENSION TO EXEMPLAR-BASED INPAINTING APPLIED TO SCRATCH CORRECTION IN DAMAGED IMAGE SEQUENCES

This paper presents an adaptation of the exemplar-based inpainting process for image sequences. The aimed application is digital restoration of old movies, and particularly scratch removal, which still remains an unsolved problem. Most of the current correction methods act frame by frame and fail to preserve the temporal coherences of the processed movies. Our new scratch concealment algorithm extends the recently published exemplar-based inpainting method by taking into account, for the block matching step, several images from the sequence. Furthermore, the method is improved by explicitely using motion estimation. A few results demonstrating the effect of this innovative extension are presented and compared to static methods.

Introduction

Achievements in the field of movie restoration are presented in this paper. Typically, movies undergo several types of degradations, caused by time and frequent usage or by accidents or other unexpected events occurred during production. Our research on movie restoration is a part of a large European project dealing with preservation and restoration of audiovisual materials 1 .

Our contribution to this project mainly concerns scratch concealment. For high-quality restoration of images or movies, a two-step method is classically used : artifacts detection (scratch or dust spots for instance), followed by artifacts removal by using specific interpolation techniques. The first step, i.e. the prior detection step for line scratches, has been published in [START_REF] Besserer | Detection and tracking scheme for line scratch removal in an image sequence[END_REF] and gives correct results. The proposed method uses projections of each image of the input sequence as input. First, a 1D-extrema detector provides candidates. Next, a MHT (Multiple Hypothesis Tracker) uses these candidates to create and keep multiple hypothesis. As the tracking goes further through the sequence, each hypothesis gains or looses evidence. To avoid a combinatorial explosion, the hypothesis tree is sequentially pruned, preserving a list of the best ones. As hypothesis are set up at each iteration, even if no information is available, a tracked path might cross gaps (missed detection or speckled scratches). Throughout this paper, we consider that we have corrupted movie data as well as scratch positions delivered by our scratch detector (so-called "mask files", see Fig. 1).

However, the scratch removal issue has not been satisfactorily solved yet in the literature. Thus, a new algorithm, belonging to the group of the inpainting techniques, is presented in this paper. It is based on an exemplar-based inpainting technique, which is extended and applied on different corrupted sequences.

Review of existing scratch correction algorithms

Scratches are very current defects in damaged film archives, characterized as bright or dark vertical lines (or colored in case of color films), spread either over several frames or being restricted in size. The temporal aspect of this degradation implies that a well performed restoration needs to be invisible in all individual frames, but also when the entire sequence is played. One major problem of the al-ready proposed correction methods lies in the fact that they are not able to reproduce the original grain of the movie when interpolating scratches, resulting in non-natural looking restored frames. This is the case of the methods based on low-pass or median filtering [START_REF] Geman | A nonlinear filter for film restoration and other problems in image restoration[END_REF][START_REF] Kokaram | Motion Picture Restoration[END_REF], for instance. Kokaram [4] used the assumption that the images follow a 2D AR model and then the sequence is generated with a 3D AR model. Despite interesting results the algorithm is, by design, not able to deal with anisotropic textures or single isophotes.

Partial success has been achieved for scratch removal in the particular case where the original texture is preserved to a certain extent. The scratches appear mostly due to abrasion by mechanical parts of the playing / recording devices and it does not always damage all the film layers. The method described in [START_REF] Joyeux | Tracking and map reconstruction of line scratches in degraded motion pictures[END_REF] uses Fourier series expansions to separate low and high frequencies of the frame data. The low frequencies (basic image structure) are reconstructed by means of a polynomial interpolation. Afterwards, a specific bayesian approach using the residual information is applied to restore the high frequencies. This algorithm performs very well and keeps the temporal coherence of the restored frames, even if the sequence is shown.

Unfortunately, the restoration is not satisfactory when texture information is entirely lost. More specific techniques belonging to the inpainting methods have to be applied. These methods allow to sequentially fill the areas of missing textures in a way that copies (or at less imitates) textures from the known areas.

Review of existing inpainting algorithms

One of the main approaches tries to prolong image contours (isophotes) from the vicinity of the missing textures [START_REF] Bertalmio | Image inpainting[END_REF]. It is related to the diffusion process in physics and PDE's (partial differential equations) are typically used. The results of these methods are convincing when removing small defects. However, diffusion PDE's cannot reconstruct textures and some blur appears on scratched and textured regions. This is particularly noticeable when dealing with old grainy films.

In [START_REF] Bertalmio | Simultaneous structure and texture image inpainting[END_REF], isophotes prolongations with PDE's and texture synthesis are combined. The image is decomposed in two sub-images, containing separately texture information and main image structure. Reconstructing separately these two images and mixing them allows to perform a textured inpainting. Nevertheless, textures in the vicinity of the isophotes still remain blurred, and macro-textures cannot be reconstructed by such a method.

Advantages of the previous approaches are combined in [START_REF] Bornard | Missing data correction in still images and image sequences[END_REF][START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF]. The missing area is filled in by an isophote-driven image sampling process. The observation and the assumption have been made that exemplar-based texture synthesis is sufficient for propagating linear isophotes within the missing data area. Very interesting results are obtained by this recent inpainting algorithm, although, so far, this algorithm has only been applied to still images. In this paper, we propose to extend it to image sequences. We investigate two ways to perform this temporal extension : the straightforward frame-by-frame application and using the estimation of motion between successive frames.

In the next section, we describe the original exemplarbased inpainting algorithm in detail and we apply it for our problem of scratch removal. Then, different temporal extensions are discussed and applied and some results are presented2 .

2 Exemplar-based inpainting applied to scratch removal Briefly, the algorithm developed in [START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF] starts the filling process at the location where the isophotes hit the boundary of the unknown area. Thus, these isophotes are prolongated first : each pixel of the boundary is characterized by a recovering priority value. The neighborhood (pattern) of the pixel with the maximum priority is then compared to all patterns of the known image area using a block-matching technique, in order to find the most similar pattern (in the sense of the L 2 -norm). Afterwards, the unknown pixels of the first pattern are replaced with the corresponding pixels from the second pattern. This process is repeated until all the unknown region is filled (Fig. 2).

Presentation of the algorithm

In this example the maximum priority is given to p and the first pattern centered on this pixel is determinated (black square). Only the known pixels of this pattern (gray part in Fig. 2) are considered for the block-matching step, which compares the first pattern to all the patterns entirely included in the known area of the image. The most similar pattern is found (white square, centered on q) and the unknown pixels of the first pattern (white area of the black square) are attributed the corresponding part of the second pattern values. Fig. 3 shows a result for such a single image (intra-frame) processing.

Modification of the exemplar-based algorithm

In order to make the algorithm more specific to our purpose of scratch removal, we propose some modifications. The algorithm presented in [START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF] can be used to fill an area of any shape. Thus, the computational load is very high, since the priority values for every pixel of the target area boundary are re-evaluated at each iteration. This is not necessary in our case : scratches often have the same geometric characteristics and the distance between points from opposite sides of the corrupted area is small, so a concentric filling of the unknown region is sufficient. All the priority values of the boundary pixels are thus calculated only once, speeding up the entire process.

We also propose another modification, related to the way the pixels are copied into the target area. In [START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF] the filling is made by blocks, resulting in the apparition of block artifacts (similar to JPEG compression). This can be avoided by copying only the central pixel of the estimated neighborhood after the pattern search.

The results presented in Figs. 4(b) and 4(d) are convincing, the algorithm performs well if the problem is to extend linear structure. Unfortunately, the results are less interesting when the isophotes exhibit some curves and discontinuities (presence on an extremity on one side of the missing data area, but not on the other, see Figs. 4(f) and4(h)). Indeed, the filling process in done by copying data from the image into the missing area. If there are no correct blocks that match the region, many visible artifacts can appear. Favorably, in our case, a sequence of images is available. In the next section, we propose to use the temporal information in order to overcome this problem.

Temporal extension

Hypothesis concerning the scratched area

A natural way to fill in the scratched area is to find correct blocks on the neighbor images when valid blocks are unavailable in the current frame. Indeed, a scratched object in a frame may be unscratched in another frame. This particularly happens when the scratch is moving and the scene is motionless or when the scene is moving and the scratch is motionless. Localization of the scratched objects in the neighbor frames should be thus the next step.

(a) (b) (c) (d) (e) (f) (g) (h)

Inpainting algorithm based on motion estimation

Let us suppose that an image at the time t contains a scratch that we would like to restore, using pixels coming from the reference images t-1 and t+1. For these images, we obtain a mask describing the position of the scratches using [START_REF] Besserer | Detection and tracking scheme for line scratch removal in an image sequence[END_REF]. The first step consists in estimating the motion between t and t-1 (backward motion), and between t and t+1 (forward motion). Afterwards the scratched pixel p with the maximum priority on the image t is found. Here are the possible cases and our proposed solutions :

1. Pixel p has an unscratched correspondent pixel either in the image t-1 or t+1. → The same block-matching as in [START_REF] Criminisi | Object removal by exemplar-based inpainting[END_REF] is performed on t-1 or t+1, but on a restricted area (the search is not exhaustive) centered on the corresponding pixel. This pattern search is performed to deal with the inaccuracy due to the motion estimation.

2. p has an unscratched correspondent pixel in both (previous, next) frames. → The search for the best block is fulfilled on the two neighbor images, and the chosen block is simply the most similar to the pattern centered in p.

3. p has no unscratched correspondent pixel. → The search is performed on the current frame, like in the original algorithm. In the next example, only two neighbor images have been used. In general, the motion between two frames is rather small, thus only a few pixels have an unscratched correspondent pixel (see Fig. 5(b)). To overcome this problem a number n of reference frames is chosen and the motion between these reference frames and the current frame is estimated. The bigger the number of the reference frames, the more scratched pixels have a correspondent pixel and the solution 1 and 2 can be applied. Our algorithm performs really well if a major part of the edges within the scratch are reconstructed using these solutions. Conversely, a high number of reference frames can also be a source of errors because of the motion. At the moment, n is chosen manually. Here is an algorithm summarizing what has been explained : while unknown pixels exist do compute the priority values of the boundary pixels for all these pixels sorted by priorities do for i = 1 to n/2 do Find a correspondent pix. in the image t -i or/and t + i if a correspondent pix. exists (unscratched) then use the solution 1 or 2 end if end for if the pix. isn't corrected yet then use the solution 3 end if end for end while Fig. 5 illustrates the filling process. As can be seen the correction is really convincing. In addition, the restoration is invisible when the sequences are played thanks to the temporal extension, as can be seen in the Figs. 6 and 7. In the first example, the scene is moving and the scratch is motionless. The other example presents the opposite case.

Conclusion

In this paper, an exemplar-based inpainting method has been applied to correct scratches in damaged movie archives. First, the method has been applied straightforwardly. Some defects in the restoration were extremely visible in some cases. Thus, the algorithm has been adapted to deal with and to take the advantage of images sequences, through a temporal extension of the previous algorithm. At the moment, results are convincing when the scene or the scratch are moving and when the scratched objects in a frame are unscratched in at least one neighbor frames. We are currently working on the other case, i.e. when the scene and the scratches are motionless within the entire shot.

Figure 1 .

 1 Figure 1. Left, a scratched image from the "lost world" sequence (1925, by courtesy of Lobster Film). Right, the result of our scratch detector superimposed over the picture. The detector outputs a binary image (mask files) showing the scratch location and extents.

Figure 2 .

 2 Figure 2. The exemplar-based inpainting process.

Figure 3 .

 3 Figure 3. Results of the algorithm in a single image. (a) shows an image with a large dirt blotch, (b) the mask which confines the unknown pixels and (c) the result.

Figure 4 .

 4 Figure 4. Scratch removal using exemplar-based inpainting.

Figure 5 .

 5 Figure 5. Interpolation of the scratch using neighbor images. (a) the scratch. (b) filling using the images t-1 and t+1, (c) images from t-2 to t+2 and (d) from t-3 to t+3. (e) the remaining pixels without inter-frame correspondents are filled using the intra-frame algorithm.

Figure 6 .

 6 Figure 6. Correction of a motionless scratch in a moving scene.

Figure 7 .

 7 Figure 7. Correction of a moving scratch in a motionless scene.

Images used for illustrating our work by courtesy of Centrimage (www.centrimage.com) and HS-ART (www.hs-art.com)

Acknowledgements J. Boldys was partially supported by the Czech Ministry of Education under the project No. 1M6798555601 (Research Center DAR).