François-Olivier Desmet

Dalil Hamroun

Gwenaëlle Collod-Béroud

Mireille Claustres

Dr Christophe Béroud
email: christophe.beroud@inserm.fr

Bioinformatics identification of splice site signals and prediction of mutation effects

Maturation of mRNA in eukaryotes is a very complex process that includes exon recognition through specific elements such as branch point motifs, 5' and 3' splice sites and splicing regulators. Mutations may affect these signals either directly by disrupting constitutive splice sites or indirectly by creating cryptic splice sites. We thus evaluated the prediction efficiency of nine software programs based on Position Weight Matrices, Markov Models, Maximal Dependence Decomposition, Neural Networks or Feature Generation Algorithms using 623 mutations for which the consequence on mRNA has been demonstrated in vitro. Position Weight Matrix-based tools were the most efficient in predicting the impact of a given mutation. Overall, at least one system correctly predicted 100% of mutations affecting invariant positions as well as -1, +3 and +5 positions of the 5'ss. Deep intronic mutations resulting in the activation of cryptic exons were almost all correctly predicted (92.31%), while other intronic mutations were less efficiently (70-80%). Exonic mutations that create cryptic splice sites were also efficiently detected (70%). Because of the prediction heterogeneity and specificity, a single tool could not be used for all predictions. Moreover, since these systems are all based on text analysis an in vitro validation step is still required.

INTRODUCTION

With the recent breakthroughs in sequencing and genotyping, millions of single nucleotide polymorphisms (SNPs) have been characterized in the human genome. Among the 7,736,157 SNPs reported in genes at NCBI (dbSNP build 130; http://www.ncbi.nlm.nih.gov/SNP/), 86.20 % are localized in introns, 2.74% in exons and the others at 5' and 3' regions [START_REF] Sherry | dbSNP: the NCBI database of genetic variation[END_REF]. Among the exonic variations, 52.19% are missense and 36.86% synonymous substitutions and may represent either neutral variations or pathogenic mutations. A review of data available from the Human Gene Mutation Database (http://www. hgmd.cf.ac.uk/ac/index.php) confirms these findings with 61.10% of diseasecausing mutations being missense mutations [START_REF] Stenson | The Human Gene Mutation Database: 2008 update[END_REF]. Concomitantly, 1610 mutations (2.55%) have also been reported that affect the invariant positions +1, +2, -1 or -2 of donor (5'ss) and acceptor (3'ss) splice sites, which delineate the limits of exons. These data are nevertheless far from reality as it is now recognized that many missense variations (up to 50% of diseasecausing single nucleotide variations) could indeed have a critical impact on mRNA maturation through disruption or creation of splice sites. It has also been shown that these variations can affect other splicing motifs such as auxiliary sequences also known as splicing regulators (Exonic Splicing Enhancers (ESE) or Exonic Splicing Silencer (ESS)) [START_REF] Zhang | RNA landscape of evolution for optimal exon and intron discrimination[END_REF][START_REF] Chasin | Searching for splicing motifs[END_REF][START_REF] Sironi | Silencer elements as possible inhibitors of pseudoexon splicing[END_REF]. In parallel, other intronic variations affect donor and acceptor splice sites although they do not concern the invariant positions. The prediction of the consequences of a nucleotide variation on mRNA is thus of major interest for molecular diagnosis not only of intronic mutations but also of exonic mutations. As the majority of disease-causing mutations that affect mRNA maturation modify donor or acceptor splices sites, many tools have been designed to predict the consequence of nucleotide substitutions on these motifs. Two types of introns have been described [START_REF] Burge | Evolutionary fates and origins of U12-type introns[END_REF][START_REF] Sharp | Classification of introns: U2-type or U12type[END_REF]. U2 snRNP-dependent introns represent more than 99.9% of all introns [START_REF] Sharp | Classification of introns: U2-type or U12type[END_REF] and are excised by a spliceosome containing the U1, U2, U4, U5 and U6 snRNPs. U12 snRNP-dependent introns are the minor class of introns and are excised by a spliceosome containing U11, U12, U4atac, U6atac and U5 snRNPs [START_REF] Boldina | Using profiles based on nucleotide hydrophobicity to define essential regions for splicing[END_REF] and are usually shorter than U2 introns. Both U2 and U12 snRNP-dependent introns are recognized by the cellular machinery through consensus motifs called donor (5'ss) and acceptor (3'ss) splice sites [START_REF] Jacob | The 5' splice site: phylogenetic evolution and variable geometry of association with U1RNA[END_REF]. These motifs are characterized by their terminal dinucleotides GT-AG, GC-AG and AT-AC, the vast majority (>98%) being of the GT-AG subtype [START_REF] Burset | Analysis of canonical and non-canonical splice sites in mammalian genomes[END_REF]. If these 4 bases are very conserved among genes, adjacent positions around the splice sites may vary, leading to short degenerated consensus sequences: usually 9 bp for the 5'ss and 14 bp for the 3'ss. Since these motifs are small, the probability to find them by chance within an intron is high thus creating pseudoexons or decoy sites. Pseudo-exons are intronic sequences that match the exon requirements, but are not selected as exons by the spliceosome. This is usually related to the splice site motifs themselves but also to the RNA structure that impacts the motif accessibility [START_REF] Buratti | Influence of RNA secondary structure on the pre-mRNA splicing process[END_REF] and/or other splicing signals. Indeed, three types of sequences are involved in splicing ( 14): constitutive splice sites, branch point se-quences [START_REF] Gao | Human branch point consensus sequence is yUnAy[END_REF] and splicing regulators (ESE and ESS), which are important for discriminating real exons from pseudo-exons [START_REF] Chasin | Searching for splicing motifs[END_REF][START_REF] Blencowe | Exonic splicing enhancers: mechanism of action, diversity and role in human genetic diseases[END_REF][START_REF] Cartegni | Listening to silence and understanding nonsense: exonic mutations that affect splicing[END_REF]. Today, research for detecting real exons and splice sites also takes into account the link between chromatin and gene structure [START_REF] Schwartz | Chromatin organization marks exon-intron structure[END_REF][START_REF] Tilgner | Nucleosome positioning as a determinant of exon recognition[END_REF] as well as hydrophobicity profiles [START_REF] Boldina | Using profiles based on nucleotide hydrophobicity to define essential regions for splicing[END_REF]. If the splicing process is complex and much work has to be done yet to fully understand it, various software and web applications are already available on Internet to detect splicing signals, at the forefront of which stands donor and acceptor splice site predictions [START_REF] Houdayer | Evaluation of in silico splice tools for decision-making in molecular diagnosis[END_REF]. To detect 5'ss and 3'ss, six methods have been developed [START_REF] Zhang | Computational searches for splicing signals[END_REF] and implemented in various software [START_REF] Brent | How does eukaryotic gene prediction work[END_REF]. The most frequently used relies on Position Weighted Matrices (PWM) and Position Specific Score Matrices (PSSM). Senapathy and Shapiro [START_REF] Senapathy | Splice junctions, branch point sites, and exons: sequence statistics, identification, and applications to genome project[END_REF][START_REF] Shapiro | RNA splice junctions of different classes of eukaryotes: sequence statistics and functional implications in gene expression[END_REF] defined PWM, which are based upon sequence alignments [START_REF] D'haeseleer | How does DNA sequence motif discovery work?[END_REF][START_REF] Zhang | A weight array method for splicing signal analysis[END_REF]. PWM have subsequently been refined and the weight of each position has been modulated to underline their respective biological importance [START_REF] Desmet | Human Splicing Finder: an online bioinformatics tool to predict splicing signals[END_REF]. Application software based on PWM processes a given sequence into short sequences (whose length is equal to that of the 5' or the 3' splice site). A weight is then attributed to each nucleotide according to its position within the PWM. If the sum of the consecutive nucleotides weights is superior or equal to a specific threshold [START_REF] Touzet | Efficient and accurate P-value computation for Position Weight Matrices[END_REF], a motif is consequently detected. In order to increase the efficiency of 5'ss detection, Carmel et al. proposed to combine a Delta-G (DG), which predicts the variation of free energy that comes from the basepairing of the 5'ss and U1 snRNP [START_REF] Carmel | Comparative analysis detects dependencies among the 5' splice-site positions[END_REF]. The Maximal Dependence Decomposition (MDD) tool relies on the assumption that splice site positions are not independent and that conditional probabilities can be determined between adjacent and nonadjacent positions. First, the algorithm calculates the amount of dependence between the consensus indicator variable and a nucleotide indicator at the position i. This value is then used to separate the dataset (constituted of aligned sequences) in two subsets: sequences with and sequences without consensus nucleotide at the given position i. This is repeated for all positions to create a tree in which each subset will be used to generate separate weight matrix models that are then combined in a composite model [START_REF] Burge | Prediction of complete gene structures in human genomic DNA[END_REF]. The third approach uses Markov models (MM). The main advantage of this technique is that it takes into account bases dependencies. A Markov model [START_REF] Eddy | Hidden Markov models[END_REF][START_REF] Eddy | Profile hidden Markov models[END_REF][START_REF] Eddy | What is a hidden Markov model?[END_REF] is defined as a suite of states (for example exon, 5' splice site, intron, 3' splice site, exon). Each state is associated to probabilities to evolve (stay in the current state or move to another one). After sequence reading, the algorithm computes the most likely suite of states (and corresponding sequences) that fulfills the various probabilities. Support Vector Machines (SVM) are classifiers that help to take a decision [START_REF] Noble | What is a support vector machine?[END_REF]. These classifiers need first to be trained using a set that contains both positive and negative data to efficiently assign the SVM parameters. Once the SVM is trained, it can be used to classify sequences. The efficiency of such algorithms is directly linked to the quality and the number of sequences used as positive and negative controls. The Feature Generation Algorithm (FGA) is an SVM derivative with an initial step to generate features. These features are groups of nucleotides (adjacent or not) that share a particular position or a range of positions. The obtained features are then used as input for a learning algorithm similar to that of SVM. Donor and acceptor sites use separate classifiers. Once the training is finished, each feature is weighted. When a sequence is processed for splice site prediction, the classifier checks if the sequence contains any feature according to the considered candidate (acceptor or donor) [START_REF] Dogan | Features generated for computational splice-site prediction correspond to functional elements[END_REF]. Neural Networks (NN) use artificial neuron networks that mimic the function of real neuron networks [START_REF] Krogh | What are artificial neural networks[END_REF]. As for SVM, a training (learning) period to set up thresholds and weights is required for the network to work efficiently. Each artificial neuron will input new pieces of information (coming from external sources or other neurons) and weight them. Once a defined threshold is reached the sequence is considered as being identified and the output will be set to 1 (0 otherwise). The last approach relies on Hybrid models. It has been shown that, in many cases, before processing data using an SVM-based approach, it could be more efficient to first sift data with Markov Models. This has been applied to splice site recognitions by Ho and Rajapakse [START_REF] Ho | Splice site detection with a higherorder markov model implemented on a neural network[END_REF], who created an hybrid approach consisting of two first-and two second-order Markov chain models followed by a threelayer neural network. Because of the strong needs to efficiently predict the effects of sequence variations on splice site recognition, not only for molecular diagnosis of human genetic diseases but also for various therapeutic approaches, a comparison of the various methods and associated algorithms/ software was necessary. In this review we evaluated various available tools/methods to predict the impact of mutations on splice site recognition using a large set of validated mutations.

MATERIAL & METHODS

Method

We mined the literature for mutations affecting splicing and selected 623 mutations, using data from UMD locus-specific databases [START_REF] Collod-Beroud | Marfan Database (third edition): new mutations and new routines for the software[END_REF][START_REF] Beroud | UMD (Universal Mutation Database): 2005 update[END_REF][START_REF] Beroud | UMD (Universal mutation database): a generic software to build and analyze locus-specific databases[END_REF][START_REF] Frederic | The FBN2 gene: new mutations, locus-specific database (Universal Mutation Database FBN2), and genotype-phenotype correlations[END_REF], the Human Gene Mutation Database (2) and other previously used datasets [START_REF] Rogan | Using information content and base frequencies to distinguish mutations from genetic polymorphisms in splice junction recognition sites[END_REF][START_REF] Rogan | Information analysis of human splice site mutations[END_REF]. Mutations were divided in four datasets (Figure 1): i) 72 mutations that affect the four invariant positions of 5' (n=49) and 3' splice sites (n=23) (positions +1, +2 and -2, -1,) (positive control); ii) 178 mutations that include intronic mutations either localized at splice sites in non-canonical positions (n=148), distant intronic mutations also known as "deep" intronic mutations (n=13) and short distance mutations (n=17); iii) 288 exonic mutations reported to affect splicing, including 10 exonic mutations that activate a cryptic splice site; and iv) 85 mutations that do not affect splice sites (negative control). These mutations have been re-ported to modify branch points (9 mutations), splicing regulators (72 mutations) or have been identified as polymorphisms in splice sites (4 mutations) during diagnosis procedures. For each set, both reference and mutant sequences were extracted from the Ensembl database version 44 (44) using Human Splicing Finder [START_REF] Desmet | Human Splicing Finder: an online bioinformatics tool to predict splicing signals[END_REF]. As some prediction programs require long intronic sequences, 200 nucleotides-long intronic sequences were added to each exon extremity. In order to facilitate data processing all sequences were stored as FASTA file (one file per set, cf. Supplementary Materials). When possible, stand-alone versions of programs were preferred because they allowed the automation of the various prediction steps. As most tools have been designed to predict the position of splice sites within a sequence mainly to detect exons, they do not allow an easy comparison of a set of mutant sequences with the corresponding reference sequences. We thus developed a Java application that combines all results in a single file and outlines the differences between sequences (reference vs. mutant) for each mutation. When thresholds had to be set, we chose to use default thresholds and parameters except for the MaxEntScan software for which the analysis was processed with two different thresholds: 0 (default threshold) and 2 as suggested by Coutinho et al. [START_REF] Eng | Nonclassical splicing mutations in the coding and noncoding regions of the ATM Gene: maximum entropy estimates of splice junction strengths[END_REF]. As the data obtained with the second threshold were more accurate (data not shown), these were subsequently used for all analyses. Note that data from MaxEntScan were extracted from the Human Splicing Finder website and not directly from MaxEntScan. Finally, when the downloaded software needed training, we used the included training sets.

Selected software programs

In order to evaluate all prediction methods, we selected the nine most widely used and freely available software programs: GenScan (30), GeneSplicer [START_REF] Pertea | GeneSplicer: a new computational method for splice site prediction[END_REF], Human Splicing Finder (HSF) [START_REF] Desmet | Human Splicing Finder: an online bioinformatics tool to predict splicing signals[END_REF], Max-EntScan [START_REF] Yeo | Maximum entropy modeling of short sequence motifs with applications to RNA splicing signals[END_REF], NNSplice [START_REF] Reese | Improved splice site detection in Genie[END_REF], SplicePort [START_REF] Dogan | SplicePort-an interactive splice-site analysis tool[END_REF], SplicePredictor [START_REF] Brendel | Gene structure prediction from consensus spliced alignment of multiple ESTs matching the same genomic locus[END_REF], SpliceView [START_REF] Rogozin | Analysis of donor splice sites in different eukaryotic organisms[END_REF] and Sroogle [START_REF] Schwartz | SROOGLE: webserver for integrative, user-friendly visualization of splicing signals[END_REF]. Three use only PWM (Table 1), two combine PWM with MM or DG, two apply a MDD approach, in one case combined with MM, one NN and the last uses an original approach called Feature Generation Algorithm (FGA).

Data management

In order to simplify interpretation of data from the three sets of mutations known to disrupt a wild type splice site, we converted all numerical information into a simple binary status: "broken" or not. The "broken" status was assigned when the software predicted that the mutation disrupted the natural splice site. This could be achieved either when the mutant splice site value was below the specific detection threshold or when the absolute variation between the reference and the mutant values was 10% higher than a specific variation that was set as previously described for HSF [START_REF] Desmet | Human Splicing Finder: an online bioinformatics tool to predict splicing signals[END_REF]. Similarly, for mutations known to create a cryptic splice site (datasets 2 and 3), we converted all numerical information into a simple binary state: "new" or not. The "new" status was assigned when the program predicted that the mutation created a splice site. This new splice site could either be near the natural splice site itself or in a deeper intronic region of the sequence, depending on the mutation localization. For dataset 4 (negative controls), both "broken" and "new" status were evaluated when a mutation was localized within a splice site or elsewhere.

RESULTS

The 72 mutations affecting the 3'ss and 5'ss invariant positions (Figure 1) were predicted to impact splice site motifs with a high accuracy (>95%) by four programs: HSF (100%), MaxEntScan (100%), SpliceView (100%), Sroogle (100%) and NNSplice (97.22%). Unexpectedly two software programs gave poor results: GenScan (45.83%) and GeneSplicer (63.89%). Overall results were more efficient for 3'ss than for 5'ss (93.72% vs. 83.90%) as reported in Table 2. Among the 148 mutations from the second dataset localized in splice site motifs, 24 affected the +3 position and 55 the +5 position, which are frequently involved in human diseases. Because these positions are degenerated and dependent on the sequence context, splice site prediction is more difficult. This was confirmed by the overall lower prediction efficiency (76.97% for the +5 position and 68.06% for the +3 position) and its heterogeneity (from 20.83% to 100% for the +3 position and from 14.55% to 100% for the +5 position) (Table 3). In addition to these two intronic key positions, position -1 (last base of the exon) also plays a major role in 5'ss recognition. Once again, the overall prediction efficiency was around 75% with high heterogeneity (from 38.10% to 100%) among tools. Only four software programs gave homogeneous results with efficiency higher than 80% for these 3 positions (HSF, MaxEntScan, NNSplice and SplicePort). For the mutations involving the other 3'ss (n=44) and 5'ss (n=25) positions, overall predictions were poorer (49.33% for 5'ss and 44.95% for 3'ss) (Table 3) and individual prediction efficiencies were even more scattered (from 4.55% to 77.27%). MaxEntScan and SplicePort gave the most accurate predictions. Then 17 intronic mutations localized at a short distance (<100 bp) and 13 deep intronic mutations (>100 bp) from the sec-ond dataset were tested. Since these mutations create cryptic splice sites, only data from the "New" status were taken into account. For the short distance mutations, results were unsatisfactory with a mean prediction efficiency of 27.45% and only Human Splicing Finder showed a detection rate higher than 70% (Table 4). With deep intronic mutations, which activate cryptic splice sites leading to cryptic exon inclusions, the efficiency increased to 60.68% and three prediction tools (Human Splicing Finder, MaxEntScan and Splice-View) displayed efficiencies higher than 80%. Concomitantly to these intronic mutations localized at a distance from constitutive splice sites, we also evaluated exonic mutations that do not directly affect constitutive splice sites but rather activate a cryptic splice site. As only few of such mutations have been studied experimentally (53-57), we collected a small set of 10 mutations belonging to this category. Human Splicing Finder, MaxEntScan and SplicePort predicted the activation of cryptic splice sites with good accuracy (70%), while the overall prediction was poor with only 44.44% accuracy (Table 5). Finally, to evaluate the proportion of exonic mutations that can affect splicing through the activation of a cryptic splice site, we selected a set of 267 mutations for which an effect on splicing has been demonstrated experimentally but without clues about the involved mechanisms. In this context, the average prediction was 23.01% with a heterogeneity ranging from 7.49% to 50.94% (Table 5). To evaluate the specificity of the various tools, we collected a set of 85 negative controls. Among them 4 polymorphisms were localized within constitutive splice sites at not highly conserved positions (-12 and +6). Five tools (GenScan, Human Slicing Finder, SplicePredictor, SpliceView and Sroogle) did not detect any impact on the corresponding splice sites as expected, while the others reported false positive predictions for 1 to 4 of these sites (Table 5). The analysis of the other 81 negative controls led to prediction of cryptic splice site activation in 17.86% of the cases with a heterogeneity ranging from 4.94% to 30.86% (Table 6).

DISCUSSION

Splice sites are key elements for exon recognition and therefore mutations leading to alteration of these signals have a strong impact on mRNA maturation and protein synthesis. It is now recognized that a wide range of mutations, which can be localized in introns but also in exons, may affect these signals either directly (disruption of constitutive splice sites) or indirectly (creation of cryptic splice sites). In order to evaluate the efficiency of prediction tools, we only evaluated the impact of mutations on splice sites discarding any effect on splicing regulators and branch points. For this purpose, we selected three sets of mutations for which the consequence on mRNA has been demonstrated in vitro and that encompass all situations: disruption of constitutive splice sites by affecting invariant or variable positions, activation of a cryptic exon by intronic mutations and exonic mutations. We anticipated that all prediction tools and methods could efficiently predict the effect of mutations involving the four invariant positions. Overall predictions were in agreement with these expectations and four software programs (Human Splicing Finder, MaxEntScan, SpliceView and Sroogle) showed an accuracy of 100%. Unexpectedly, GeneSplicer and Genscan gave poor results especially for 5'ss. Three other positions play an important role in 5'ss recognition: the -1, +3 and +5 positions. It has been recently demonstrated that the +3 position is intrinsically linked to sequence contexts leading to difficult predictions [START_REF] Guedard-Mereuze | Sequence contexts that determine the pathogenicity of base substitutions at position +3 of donor splice-sites[END_REF], while the +5 position is less influenced [START_REF] Krawczak | Single base-pair substitutions in exon-intron junctions of human genes: nature, distribution, and consequences for mRNA splicing[END_REF]. Interestingly, despite these limitations MaxEntScan could accurately predict all pathogenic mutations affecting these positions. HSF, NNSplice and SplicePort tools also gave very accurate results. Like for the "invariant" positions, GenScan gave poor predictions while GeneSplicer, SplicePredictor and SpliceView showed intermediate accuracy. Sroogle was very efficient for the -1 and +5 position (100%), but performed less efficiently for the +3 position (54.17%). This is probably due to the use of a PWM that does not take into account the sequence context, this limit not being compensated by the addition of the DG algorithm. Besides these mutations that disrupt a constitutive splice site, other mutations can activate a cryptic splice site that will either be in competition with the constitutive site for recognition by the cellular machinery or lead to the inclusion of a pseudo exon [START_REF] Tuffery-Giraud | Pseudoexon activation in the DMD gene as a novel mechanism for Becker muscular dystrophy[END_REF][START_REF] Beroud | Dystrophinopathy caused by mid-intronic substitutions activating cryptic exons in the DMD gene[END_REF]. Since these sites are localized in introns at short or long distance from the constitutive splice sites, they are usually harder to predict, as pseudo exons do not match criteria for exon recognition by Markov Models and Neural Networks. On the other hand, PWM-based tools should accurately identify them because they analyze only splice sites motifs. As expected PWM tools gave the best results especially for the deep intronic mutations (average of 73.85%), but performed less well in the identification of intronic mutations localized at a short distance from constitutive splice sites (average 30.59%). Only HSF gave accurate results in both situations (92.31% and 70.59%). Programs based on Markov Models and Neural Networks were less efficient with a detection rate ranging from 0 to 69.23% with poor results for short distance cryptic sites (5.88 to 41.18%). Exonic mutations can also result in cryptic splice site activation. Because these sites are localized in a favorable context for recognition by Markov Models and Neural Networks, we expected them to be as efficient as PWM. Indeed, the best NN-and Table 6: Prediction of intronic and exonic mutations that do not affect splicing. ss polymorphism = mutations localized in constitutive splice sites; other positions = intronic or exonic mutations not localized in constitutive splice sites.

PWM-based tools gave similar results with a detection rate of 70% (HSF, Max-EntScan and SplicePort) while the other three PWM-based tools gave unexpectedly low predictions, especially Sroogle (0%). MDD-based tools were less efficient as for previous situations. In order to evaluate the proportion of exonic mutations that affect splicing and are associated with splice site modifications (disruption of the constitutive splice site or creation of a cryptic splice site) we selected a series of 267 exonic mutations belonging to this category. If we consider that only 70% of exonic mutations that result in the activation of a cryptic exon are correctly predicted by the most efficient tools (HSF, MaxEntScan and SplicePort) and we remove tools for which a false positive prediction rate is high (SplicePort, MaxEntScan) we can predict that between 31 to 44% of exonic mutations may alter splicing by involving a splice site. Finally, to evaluate the specificity of predictions we selected 85 intronic and exonic mutations for which it has been demonstrated that they do not affect constitutive or cryptic splice sites (negative control). Four of them were localized in constitutive splice sites and therefore directly reflected false positive predictions for the inactivation of these sites. Five tools (GenScan, HSF, SplicePredictor, Splice-View and Sroogle) did not predict any consequence for these mutations. Conversely, NNSplice predicted the disruption of one of the constitutive splices sites, GeneSplicer and MaxEntScan of two and SplicePort of all. If these results are not significant due to the small number of mutations, they question the specificity of these last three tools. The evaluation of predictions for the other exonic and intronic mutations revealed a low rate (17.86%) of positive predictions. Because the creation of a cryptic splice site is not sufficient to activate a cryptic exon, these results could not be considered as false positive. In fact, two main situations could be encountered. The cryptic splice site is in the vicinity of the constitutive splice site and therefore is in competition for recognition by the cellular machinery. The strength of this new site, its localization in relation to the branch point for 3'ss as well as the sequence context (splice regulators) should all be taken into account before predicting if this splice site will be recognized by the cellular machinery. To our knowledge, today no tool can handle such information and a manual analysis has to be performed. In the second situation, the cryptic splice site is at a distance from an exon. In this context a cryptic exon can be recognized only if the complementary splice site (3'ss for a 5' cryptic ss and vice versa) is present at a short distance. In addition, other criteria should also be satisfied: presence of a branch point and a favorable splice regulator context. Once again, the tools tested in this study can not handle such information. GeneSplicer and GenScan were expected to perform better but cryptic exons are usually different from natural exons (no selection for codons and presence of stop codons), reducing their ability to recognize the cryptic exons.

Predictions limitations

Several limitations should be taken into account when interpreting the predictions of the tools used for this analysis. The first is related to the intron type. Although two types of introns have been described [START_REF] Burge | Evolutionary fates and origins of U12-type introns[END_REF][START_REF] Sharp | Classification of introns: U2-type or U12type[END_REF], only the GT-AG subtype of U2 snRNP-dependent intron model is used for predictions as it accounts for more than 98% of the cases. Consequently, a specific attention should be paid to wild type splice sites to confirm that they belong to this category before performing prediction analyses. The second limit is associated with the algorithms. They are all based on a text analysis (search or comparison) without addition of other parameters such as accessibility of the motif based on 2D or 3D structures as well as DNA-protein interactions that could mask a site. Thus strong splice sites could be predicted but not used in vivo because they are not accessible to the spliceosome. Moreover, the splicing process is complex and involves many signals (branch point, splicing signals, splice regulators) and proteins, which are not taken into account for predictions [START_REF] Nilsen | The spliceosome: the most complex macromolecular machine in the cell?[END_REF]. The software programs we tested in this study use different methods to achieve splice site recognition. They all rely on the quality of their respective training/building sets. With the completion of the Human Genome Project, many exons have been identified and are available to build positive sets of 5' and 3'ss. This is particularly useful for PWM-based tools that only require a positive control set. On the contrary, for SVM methods, both positive and negative control sets are required. In addition, some software programs need to know where the exon-intron boundaries are while others make ab initio prediction. In this review, we tested both splice site detectors and gene/exon annotation programs. The best results were obtained with PWMbased software. Gene/exon annotation programs (GenScan, GeneSplicer) have a global approach and need a longer sequence to detect splice sites (ranging from 80 to 200 intronic nucleotides) and evaluate splice regulators to efficiently detect true exons. This additional step is critical and thus led to a reduced efficiency of such tools in this analysis. For instance, GenScan had first to detect exon boundaries, which was not possible for all reference sequences. This does not impair their ability to predict exons in a largescale genomic analysis.

Which software for which situation?

To decide which is the most appropriate software, the user has to consider not only its efficiency, but also additional characteristics such as its availability, amount of data to be processed (batch analysis or multiple query option), interface, etc. In this review we have tested the nine most frequently used tools. Only Human Splicing Finder was created to compare a mutant sequence with a reference sequence, while the other tools were designed to identify splice sites within a sequence (MaxEntScan, Sroogle, SplicePort, SplicePredictor and SpliceView) or to de-tect exons (GeneSplicer and GenScan). Because of these different designs and underlying algorithms, they displayed different efficiencies. To evaluate a mutation involving invariant positions most tools may be used with a preference for Human Splicing Finder, MaxEntScan, NNSplice, SpliceView and Sroogle, the last two being less efficient when positions -1, +3 and +5 of 5'ss need to be evaluated. When other less conserved positions have to be investigated, users may prefer Max-EntScan, NNSplice or SplicePort. Nevertheless, it is also critical to evaluate the specificity of the different prediction tools as MaxEntScan, NNSplice or SplicePort seem to have a reduced specificity compared to Human Splicing Finder, SplicePredictor, SpliceView and Sroogle. When searching for a cryptic splice site, Human Splicing Finder is the most efficient, but MaxEntScan and SpliceView also are very efficient for deep intronic mutations. Finally, for exonic mutations Human Splicing Finder, MaxEntScan, SplicePort and NNSplice are the most efficient. In this study, we focused on splice site detection, one of the many features of the complex mechanism of splicing. In diagnostic and research situations, geneticists have to address all aspects of splicing defects including alteration of branch points and the effect of splicing regulators since they modulate splicing by allowing fixation of proteins/complexes on the pre-mRNA, thus enhancing weak splice site signals or decreasing strong splice site signals. The user should thus combine the tools presented here with other systems/approaches to evaluate these many situations [START_REF] Zhang | RNA landscape of evolution for optimal exon and intron discrimination[END_REF][START_REF] Sironi | Silencer elements as possible inhibitors of pseudoexon splicing[END_REF][START_REF] Desmet | Human Splicing Finder: an online bioinformatics tool to predict splicing signals[END_REF][START_REF] Cartegni | ESEfinder: A web resource to identify exonic splicing enhancers[END_REF][START_REF] Goren | Comparative analysis identifies exonic splicing regulatory sequences--The complex definition of enhancers and silencers[END_REF][START_REF] Zhang | Computational definition of sequence motifs governing constitutive exon splicing[END_REF][START_REF] Wang | Systematic identification and analysis of exonic splicing silencers[END_REF][START_REF] Fairbrother | RESCUE-ESE identifies candidate exonic splicing enhancers in vertebrate exons[END_REF]. Although Human Splicing Finder is not the most efficient tool for all situations, it is the first that can address all splicing aspects (branch points, splice sites and splicing regulators). Users can thus include in their toolbox Human Splicing Finder, MaxEntScan and Sroogle. Finally, since all in silico predictions are based only on statistics, an in vitro validation step is still required.

Figure 1 :

 1 Figure 1: Distribution of mutation sets according to splicing motifs: 5'ss (donor splice site); 3'ss (acceptor splice site); branch point and splicing regulators (ESE, ESS).

Table 1 :Table 2 :

 12 Selected software programs to predict splice sites. PWM = Position Weight Matrix; MDD = Maximal Dependence Decomposition; MM = Markov Model; NN = Neural Network; FGA = Feature Generation Algorithm; DG = Delta-G Prediction of mutations affecting the invariant positions -2, -1, +1 and +2 of 3'ss (n=23) and 5'ss (n=49).

Table 3 :

 3 Prediction of mutations affecting the last base of the exons as well as positions +3 and +5 of the 5'ss and of mutations affecting other positions of 3'ss and 5'ss.

Table 4 :

 4 Prediction of intronic mutations localized at a distance from splice sites.

Table 5 :

 5 Exonic mutations known to result in splice defects.

FUNDING


Charalampos Karypidis
email: ch_karypidis@yahoo.com

ORDER EFFECTS AND VOWEL DECAY IN SHORT-TERM MEMORY: THE NEUTRALIZATION HYPOTHESIS

Keywords: order effect, vowel perception, neutralization, short-term memory, peripherality

Information stored in short-term memory decays extremely fast compared to the content of long-term memory. The nature of this memory loss being little known (random or systematic), this paper examines the presentation order effect in the light of the neutralization hypothesis, according to which the first vowel in a pair decays, while retained in short-term memory, toward [@]. Twelve French listeners participated in three AX roving discrimination sessions. For each of the 10 vowel categories of French, a prototype and four satellite tokens were synthesized. Tokens were paired with themselves as well as with the prototype (in both orders of presentation). Results revealed minor or major order effects within every phonetic category. However, the neutralization hypothesis could not account for at least half of these asymmetries. An alternative model, which adopts the periphery of the vowel space as a reference area, appears to fit best with our data.

INTRODUCTION

The 'presentation order effect' (POE) -also known as 'phenomenon of asymmetry' -can be summarized in Polka and Bohn's [START_REF] Polka | Asymmetries in vowel perception[END_REF] definition: "Asymmetries in vowel perception occur such that discrimination of a vowel change presented in one direction is easier compared to the same change presented in the reverse direction".

The order effect has been attributed to various factors such as peripherality (see Polka [START_REF] Polka | Asymmetries in vowel perception[END_REF] for F1/F2; Best [START_REF] Best | Developmental increase in infants' discrimination of nonnative vowels that adults assimilate to a single native vowel[END_REF] for F1/F2/F3), focalization [START_REF] Schwartz | Asymmetries in vowel perception, in the context of the Dispersion-Focalisation Theory[END_REF] and typicality [START_REF] Kuhl | Human adults and human infants show a "perceptual magnet effect" for the prototypes of speech categories, monkeys do not[END_REF]. Using Rosch's (8) terminology, we will consider the aforementioned factors as cognitive reference points.

Cowan and Morse (3), while interpreting their results, suggested that in a pair of vowels V1V2, it is possible that the auditory trace left in memory by V1 (the first vowel presented in the pair) might gradually degrade toward a neutral (central) position in the vowel space which acts as a perceptual anchor. The neutral vowel would thus behave as a perceptual magnet toward which a vowel is drawn whilst stored in short-term memory, waiting to be used or to be discarded in order to make room for the next chunk of information to be stacked away.

Figure 1 depicts the direction of decay of stimuli 1 and 2 presented in both orders (1-2 and 2-1). According to the neutralization hypothesis, when the tokens are presented in the order 1-2, stimulus 1 will degrade toward [@] and, at the same time, toward 2, which is located on the axis of decay. Therefore, the initial acoustic difference will become smaller. In the order 2-1, on the other hand, stimulus 2 will move toward [@] but this time further away from 1, rendering the initial difference greater. Therefore, discrimination scores in the order 2-1 will be greater than those for the order 1-2. More generally, discrimination would be easier when in a pair, the first stimulus presented is closer to the adopted reference point. Repp and Crowder (7), after constructing an innovative set of stimuli to examine this hypothesis (which they called the 'neutralization hypothesis', henceforward NH), came to the conclusion that vowels tend rather to be drawn toward the center of a given stimulus set.

This paper examines the NH for French vowels. A set of 50 stimuli spanning the ten vowel categories in French allowed us to observe the direction in which vowels are easier to discriminate.

Stimuli

Nine monophthongal vowels, [i, e, E, u, o, O, a, y, ø], corresponding to the average values 1 of French vowels uttered by adult male speakers, were synthesized. The neutral vowel [@] was then added to the original set of 9 vowels. Its properties were assumed to be identical to those of a uniform tube: F 1 = 500, F 2 = 1500, F 3 = 2500, F 4 = 3500 Hz. Raw values for all ten vowels were converted into Bark scale [START_REF] Traunmüller | Analytical expressions for the tonotopic sensory scale[END_REF] in order to account for discrepancies in sensitivity of frequency perception (Table 1). Each of the ten original prototypes 2 (P) was surrounded by four neighboring stimuli (N1-N4) in the form of a cross. N1-N4 were positioned on the endpoints of each cross, one axis of which pointed toward [@]. Each arm was equivalent to a Euclidean distance (on the F1-F2 plane) of 0.4 Bark between the prototype and each neighbor. N1 was located on the axis pointing toward [@] and was positioned the furthest away from it (contrary to N3 which was positionned the closest to [@]). The rest of the tokens were numbered in a clockwise fashion. N1-N4 for the /@/ category were arbitrarily numbered and the two axes were parallel to the F1-F2 coordinates. An F1-F2 plot of the fifty stimuli is found in Figure 2.

F3 and F4 were fixed independently within each category (the F3 and F4 values of each prototype were used for its four corresponding satellites). Duration of tokens was fixed at 250 ms and F0 contour was falling (100-80 Hz). A cascade formant synthesizer (4) was used for the preparation of the stimuli, which were matched in RMS energy at -10 dB using Sound Forge 6.0.

Procedure

For each phonetic category, tokens (P, N1-N4) were paired with themselves (P/P, N1/N1, . . . ) and the four neighbors were also paired with the prototype in both orders (P/N1, N1/P, . . . ). The Inter-Stimulus Interval (ISI) was fixed at 500 ms. Six experimental blocks, each containing all pairs in random order (identical for all listeners), were prepared. The first block was considered a training session.

Listeners were requested to judge whether the paired stimuli were absolutely identical or even slightly different by typing 'd' for 'différents' (eng. different) or 'm' for 'mêmes' (eng. same) on the keyboard. The first stimulus of the following pair was presented with a 1000-millisecond delay. No feedback was given after each answer.

Results

According to NH, four main predictions (Pr) could be made (Figure 3):

• Pr1: positive order effect for N1 (P/N1>N1/P): discrimination scores for pairs in the order 'P/N1' would be greater than scores in the order 'N1/P'; • Pr2: negative order effect for N3 (P/N3<N3/P); • Pr3: no particular order effect for N2 and N4; • Pr4: no particular order effect for the [@] category. Two separate three-way (factors: Order, Neighbor, Vowel) ANOVAs were conducted, one on pairs involving N1 and N3 and another for N2 and N4 pairs. Scores for the [@] category were excluded, because no order effects were expected.

In interactions. The effect of Order for N2 and N4 was evidently incompatible with NH. At the same time, the ANOVA on N1 and N3 could not inform us whether the direction of order effects was the one predicted by NH (cf. Figure 3). Therefore, a follow-up analysis (separate ANOVAs for each Neighbor) was conducted.

This time, a significant effect of Order was found for N1 [F [START_REF] Best | Developmental increase in infants' discrimination of nonnative vowels that adults assimilate to a single native vowel[END_REF]198)

Discussion

A better understanding of the asymmetries found in our data can be obtained with a graphical representation. In Figure 4, the arrows depict the order in which listeners discriminated considerably better.

We have chosen as a threshold of significant difference between two orders that proposed by Repp and Crowder, that is, 10%. We have later readjusted this threshold to 11%, given that some scores were critical (10.8%). Therefore, for a given pair (i.e. P/N1), an arrow pointing toward N1 suggests that discrimination was better in the order 'P/N1': (P/N1)-(N1/P)> 11%]. As it was explained in the Introduction section, this is due to the fact that in this order, P is attracted toward [@] and thus distances itself even more from N1.

NH predicts 18 order effects (Figure 3), two (N1 and N3) for each of the nine vowel categories of which only 11 are found in Figure 4. In addition, NH cannot account for the 14 additional asymmetries involving N2, N3 or N4. It is also worth noting that within the /@/ category, order effects were found for 3 of four neighbors (Pr4). Practically all arrows exhibit a rather considerable consistency: almost all of them point toward the periphery (the edges) of the vowel space. This seems to agree with a series of papers on the role of peripherality in the asymmetry effect (see [START_REF] Polka | Asymmetries in vowel perception[END_REF] for a review).

It appears thus that discrimination is easier when a peripheral vowel is presented second. In this case, one is to assume that the reference point is not a single point on the space but rather involves the whole perimeter (periphery) of the vowel space. However, if the asymmetry effect was due to the vowel decaying towards a reference point (the periphery), discrimination would be easier when the vowel closer to the reference point was presented first (see Figures 1 and 3). Arrows in Figure 4 exhibit an inversed polarity. This led us to the hypothesis that the presentation order effect is rather triggered by a constrast effect due to recency effects. Our postulate would thus be: "the more peripheral the second vowel, the stronger the constrast it generates with the first vowel".

Our data and our hypothesis on the role of peripherality are supported by Repp and Crowder's results, where the role of peripherality appears to account for more than 85% of the order effects. Moreover, the direction of the arrows found in Figure 4 is very close to that presented in Polka and Bohn's (6) Figure 1a, which offers an overview of the literature on infant data.

CONCLUSION

In this paper, we have examined the presentation order effect in the light of the neutralization hypothesis. Using a set of stimuli spanning the ten phonetic categories of French vowels and three AX roving discrimination tests, we have studied the extent to which the neutralization hypothesis can account for presentation order effects obtained throughout the entire vowel space. Results indicated that the aforementioned hypothesis does not fit well with our data. A graphical representation of the order effects revealed that peripherality appears to be a much more plausible predictor. A hypothesis was then stated that asymmetries are not due to vowel decay but that are rather triggered by the stimulus presented second in a pair. More precisely, in a pair of front vowels, the order effect would be triggered by the vowel that is more front (higher F2) via a retroactive contrast effect. The same tendency would be valid for high, low and back vowels.

Whether peripherality is a factor setting off asymmetries, a very important, two-fold question remains to be answered: which cognitive process is behind this phenomenon and whether the trigger is the first or the second vowel. NH suggests that the decisive factor is the direction of decay of V1 whilst the hypothesis of peripherality presupposes a retroactive contrast effect. Any attempt to answer these questions could have important theoretical implications in the field of speech perception, giving us an insight to how information is analyzed and stored in short-term memory.

Further experiments to address this question are currently being conducted. Only pairs containing the prototype and N1 of each category were included. In the first two tasks, the Inter-Stimulus Interval was equal to either 200 or 1000 millisecond. At a second stage, the three tasks described in this paper were undertaken by listeners with different linguistic background. Given that the periphery of the vowel space is, due to articulatory and acoustic contraints, similar to all listeners, order effects from both language groups (French and non-French) are expected. For lack of space, the effect of experimental block and reaction times collected during our experiments will be exposed in a future paper.

Figure 1 :

 1 Figure 1: Direction of decay predicted by the neutralization hypothesis for stimuli 1 and 2 in the order 1-2 (a) and in the order 2-1 (b). In (a), 1 moves toward [@] and, at the same time, toward stimulus 2; in (b), stimulus 2 moves toward [@] but this time, further away from 1.

Figure 2 :

 2 Figure 2: F1-F2 plot of the fifty synthesized stimuli. In the center of each cross lies the prototype. Neighbors (N) 1 and 3 lie on the axis pointing toward [@].

Figure 3 :

 3 Figure 3: Order effects predicted by the neutralization hypothesis. Discrimination is easier when the stimulus presented first is located closer to [@].

  = 44.325, p < .001], N2 [F (1, 198) = 9.936, p < .005] and N4 [F (1, 198) = 9.627, p < .005] but not for N3 [F (1, 198) = .186, p = .667]. The effect for N1 confirms Pr1. However, the order effect for N2 and N4, on one hand, and the absence of effect for N3, on the other, are incompatible with Pr3 and Pr2 respectively. The effect of Vowel was significant in all cases [F (8, 198) = 10.888; F (8, 198) = 7.071; F (8, 198) = 5.783; F (8, 198) = 13.259 for N1, N2, N3 and N4 respectively; p < .001 for all four neighbors].

Figure 4 :

 4 Figure 4: Order effects for the 12 subjects. Arrows point to the direction in which vowel change was perceived more efficiently.

Table 1 :

 1 Formant values of the ten vowel prototypes in Bark.

	Vowel F1	F2	F3	F4
	i	3.11 13.22 15.63 16.49
	e	3.70 12.90 14.90 16.40
	E	5.20 12.00 14.60 16.30
	a	6.40 9.90 14.50 16.20
	O	5.20 8.50 14.20 16.20
	o	3.90 7.20 13.90 16.20
	u	3.20 7.00 13.10 15.90
	y	3.00 12.10 13.40 16.00
	ø	3.80 10.70 13.80 16.10
	@	4.90 11.10 14.50 16.70
		2. EXPERIMENT
	2.1. Participants
	In total, 18 French-speaking listeners (range: 21-47
	years; mean: 28 years) participated in this experi-
	ment. Not all participants underwent all three tasks,
	thus 12 listeners were finally obtained for each of
	the three session. All reported being native speakers
	and having no known speech or hearing disorders.

In Calliope[START_REF] Calliope | La parole et son traitement automatique[END_REF], [e, o, u, y, ø] were uttered in a pV context and [i, E, a, O] in a pVö context. Nonetheless, it is the only reference, of which we are aware, that comprises F3 and F4 values. Formant values obtained from preliminary recordings of French vowels in isolated context concur with the aforementioned values.

[START_REF] Calliope | La parole et son traitement automatique[END_REF] We use the term 'prototype' here only for practical reasons. We do not presume that all listeners necessarily share the same category exemplars.

ACKNOWLEDGMENTS

We would like to thank our listeners for their time and energy as well as . . .


Trang C Mai

Malcolm Egan

Trung Q Duong

Marco Di Renzo

Event Detection in Molecular Communication Networks with Anomalous Diffusion

A key problem in nanomachine networks is how information from sensors is to be transmitted to a fusion center. In this paper, we propose a molecular communicationbased event detection network. In particular, we develop a detection framework that can cope with scenarios where the molecules propagate according to anomalous diffusion instead of the conventional Brownian motion. We propose an algorithm for optimizing the network throughput by exploiting tools from reinforcement learning. Our algorithms are evaluated with the aid of numerical simulations, which demonstrate the tradeoffs between performance and complexity.

I. INTRODUCTION

Networks consisting of a large number of nanomachines, which are able to sense, communicate, and actuate at the nanoscale, have been proposed for applications ranging from intrabody health monitoring to pollution control [START_REF] Nakano | Molecular communication and networking: opportunities and challenges[END_REF]. A key component of these networks is the detection of events such as the presence of undesirable chemicals in the atmosphere or the malfunction of cells in biological systems. In a network of distributed nanomachines, this information must then be sent to a fusion nanomachine (FN), which can take action to mitigate the effect of the event.

Molecular communication forms one approach to support communication among sensing nanomachines (SNs) and a FN, where unlike conventional electromagnetic-based communication, each SN encodes information in the release time, number, or type of molecules emitted by each SN [START_REF] Nakano | Molecular communication and networking: opportunities and challenges[END_REF]- [START_REF] Farsad | A comprehensive survey of recent advancements in molecular communication[END_REF]. Molecular communication raises new challenges, due to the low energy and limited computational resources available at each nanomachine. Moreover, the achievable throughput of molecular communication systems is limited by the noise introduced from the random diffusion time of each molecule.

In existing works, the diffusion is typically modeled according to Brownian motion. In this case, it is possible to derive closed-form expressions for the first passage time distribution and the impulse response for reasonable boundary conditions. Under the assumption of Brownian motion, capacity characterizations [START_REF] Srinivas | Molecular communication in fluid media: The additive inverse gaussian noise channel[END_REF], [START_REF] Egan | Varianceconstrained capacity of the molecular timing channel with synchronization error[END_REF], practical receiver designs [START_REF] Chiu | Near-optimal low complexity receiver design for diffusion-based molecular communication[END_REF], [START_REF] Li | convexity inspired low-complexity noncoherent signal detector for nanoscale molecular communications[END_REF] and intersymbol interference mitigation strategies [START_REF] Noel | Improving receiver performance of diffusive molecular communication with enzymes[END_REF] have been proposed.

However, Brownian motion cannot capture anomalous diffusion, where the growth of the mean square displacement is much faster or slower as it usually occurs in scenarios with turbulence [START_REF] Klafter | Beyond Brownian motion[END_REF]. A key class of models for anomalous diffusion are the Lévy walks, which introduce challenges for the design of reliable molecular communication systems, as the first passage time and the impulse response are not generally available in closed-form. Moreover, the diffusion process cannot be encapsulated in a small number of parameters in Lévy walk models, as is the case in Brownian motion. These challenges are exacerbated in the context of event detection, where the large number of SNs lead to a corresponding number of diffusion channels, each with a different distance and fluid properties between the transmitting and receiving nanomachines. In particular, existing molecular communication techniques for Brownian motion and a single transmitting SN cannot be directly applied.

In this paper, we propose a molecular communication system for event detection which can cope with anomalous diffusion. Our system consists of a event detector and a training phase to optimize the time slot duration and tune the detection rule. In particular, we design the training phase to estimate the parameters required for reliable detection based on one bit messages from each SN. To design the training phase, we introduce a reinforcement learning-based algorithm for optimizing the time slot duration. In the proposed algorithm, the duration of each time slot is chosen in order to maximize the throughput and is applicable to any diffusion process with stationary increments and an independent first passage time for each link between a SN and the FN. In particular, it can be applied to the Lévy walks. To evaluate the performance of the system, we derive a tractable expression for the probability of error and numerically investigate tradeoffs between performance and complexity.

The remainder of this paper is organized as follows. In Section II, the system model is described. In Section III, the detection rule and the error probability are introduced. In Section IV, the training phase design is provided. In Section V, the numerical results are given. Finally, Section VI concludes this paper.

II. SYSTEM MODEL

We consider a synchronized and time-slotted nanoscale network, where N S SNs monitor the state of a system. The duration of each time slot is T S . When a SN observes an event, it emits a molecule at the beginning of the next transmission slot. The intended receiver of the message is the FN, which processes the messages from the SNs. For example, the FN might be the head of a cluster of SNs which sends the detection of an event to an actuator nanomachine.

Events occur with probability p E and we assume that each SN detects an event with probability p D . The molecule transmitted by each SN diffuses to the FN according to a Lévy walk [START_REF] Zaburdaev | Lévy walks[END_REF]. We assume that the FN perfectly absorbs any detected molecule. As we will see in further text, our system design is applicable to any network of SNs where the diffusion process is stationary and the first passage time of each SN-FN link is independent. In particular, the proposed design applies both to Brownian motion and to Lévy walks under a wide range of boundary conditions. As Lévy walks form a very general class of diffusion processes, to make our description concrete, we assume that the FN is an absorbing boundary, the Lévy walk is one-dimensional, and the distance that a molecule travels within time t is given by

d t = nt j=1 vD j , (1)
where n t is the number of jumps within time t, D j ∈ R is the duration and direction of the j-th jump, and v is the fixed velocity. We assume that D j is an α-stable random variable with exponent 1 < α < 2, scale parameter σ ∈ R >0 , skewness β = 0, and shift δ = 0.

Molecules that are emitted in earlier time slots and that are not absorbed by the FN persist. As a consequence, there is inter-symbol interference and the number of molecules observed by the FN in time slot j is given by

j k=1 Y j k , where Y j k = N S i=1 Y j k,i
is the number of molecules that are transmitted in time slot k and that arrive in time slot j. In particular, the random variable Y j k,i is Bernoulli distributed with success occurring when an event occurs and the molecule sent from the i-th SN in slot k arrives at the FN in slot j. Let F (t) be probability that a molecule transmitted at time 0 is absorbed before time t. Then the success probability is equal to

P j-k,i = F ((j + 1 -k)T S ) -F ((j -k)T S ), (2)
where P j-k,i is the probability that (j -k) time slots for the molecule from SN i to diffuse to the FN are needed. Based on these assumptions, the inter-symbol interference is given by j-1 k=1 Y j k and the desired signal is Y j j .

III. DETECTION RULE AND ERROR PROBABILITY

To determine whether or not an event has occurred, the FN solves a hypothesis-testing problem based on the number of observed molecules. Let H 1 be the hypothesis that an event occur and H 0 be the hypothesis that an event did not occur. The number of molecules that arrive in time slot j under each hypothesis is given by

H 0 : R H0 = j-1 k=1 Y j k H 1 : R H1 = Y j j + j-1 k=1 Y j k , (3)
For an observation of r molecules, the detection rule is given by the likelihood ratio test [START_REF] Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF] Pr(r|H 1 ) Pr(r|H 0 )

H1 H0 η, (4)
where η = (1 -p E )/p E and the costs are C 00 = C 11 = 0 and

C 01 = C 10 = 1.
As observed in [START_REF] Chiu | Near-optimal low complexity receiver design for diffusion-based molecular communication[END_REF], Y j k,i , which is the event that the molecule emitted in time slot k by SN i that arrive in time slot j, is a Bernoulli random variable with success probability P j-k,i . It then follows that Y j k = N S i=1 Y j k,i is binomial distributed. For finite j, R H0 and R H1 , which are sums of independent binomial distributions, are Poisson-binomial distributed. Moreover, for systems operating for long periods of time corresponding to j → ∞, it follows from the law of rare events [START_REF] Chen | On the convergence of Poisson binomial to Poisson distributions[END_REF], [START_REF] Cam | An approximation theorem for the Poisson binomial distribution[END_REF] that R H0 and R H1 converge in distribution to a Poisson random variable. The large j approximation also forms a worst case scenario for short system run times.

To derive the parameter of the Poisson random variable, we note that an event occurs and is detected by SN i with probability p E p D . Moreover,

j k=0 P j-k,i = 1 when j → ∞.
This means that the inter-symbol interference due to SN i is a Poisson random variable with expected value p E p D (1 -P 0,i ), where P 0,i is the probability that a single time slot is needed for a molecule emitted by SN i to reach the FN. Since the sum of Poisson random variables is Poisson, it follows that

R H0 ∼ P oisson p E p D N S i=1 (1 -P 0,i ) R H1 ∼ P oisson N S i=1 P 0,i + p E p D N S i=1
(1 -P 0,i ) . [START_REF] Egan | Varianceconstrained capacity of the molecular timing channel with synchronization error[END_REF] We are now in the position to derive a decision rule. Let

λ 0 = p E p D N S i=1 (1 -P 0,i ) λ 1 = N S i=1 P 0,i + p E p D N S i=1 (1 -P 0,i ). (6)
By using the probability density function of the Poisson distribution in (4), we obtain the decision rule

r H1 H0 Θ, (7)
where

Θ = log η + λ 1 -λ 0 log(λ 1 ) -log(λ 0 ) . ( 8
)
Under the Poisson approximation, the average probability of detection error is given by

P e = p E F (λ 1 ; Θ ) + (1 -p E )(1 -F (λ 0 ; Θ )), (9)
where F (λ; x) is the CDF of the Poisson distribution. Since Poisson-binomial random variables converges in distribution to a Poisson random variable under the conditions [START_REF] Chen | On the convergence of Poisson binomial to Poisson distributions[END_REF] satisfied by our model, it follows that the probability of error converges to [START_REF] Klafter | Beyond Brownian motion[END_REF].

A key observation is that the decision rule and hence the error probability are completely determined by the parameter

A = N S i=1 P 0,i . (10)
We remark that if Brownian motion is considered, it is possible to obtain an analytical expression for P 0,i , under appropriate boundary conditions, and perfect knowledge of the fluid diffusion coefficient as well as the distances between each SN and the FN. However, in the case of Lévy walks, a closed-form expression for the first passage time is not available. Moreover, due to the large number of SNs, it is unlikely that the FN will have access to the necessary information on the fluid properties and locations of all the SNs throughout the network. In the next section, we explore the effect of the training phase on the estimate of A and the decision rule.

IV. TRAINING PHASE DESIGN

Due to the unknown diffusion process, the training phase plays a crucial role in obtaining the parameters that are required for the decision rule. In this section, we design the training phase to ensure that the probability of error lies within a known confidence interval. We then exploit the confidence interval to develop an online algorithm inspired by reinforcement learning [START_REF] Kaelbling | Reinforcement learning: A survey[END_REF], to optimize the time slot duration. By exploiting reinforcement learning, our approach optimizes the time slot duration by exploiting observations from sensor signals without having a complete characterization of the channel. As the first passage time and the impulse response in our model are not available in closed-form, such a learningbased approach is highly desirable.

To develop an online algorithm, we assume that after each time slot in the training phase, the network can remove the transmitted molecules from the system. This can be achieved by the introduction of reactants, such as the enzyme-based approach in [START_REF] Noel | Improving receiver performance of diffusive molecular communication with enzymes[END_REF], which reduce the lifetime of the information molecules, only necessary during the training phase to reduce the energy consumption at the FN while achieving accurate parameter estimation.

From [START_REF] Zaburdaev | Lévy walks[END_REF], it follows that the probability of error is completely determined by the parameter A in [START_REF] Zaburdaev | Lévy walks[END_REF]. As a consequence, the key step in the training phase is to estimate A. A natural estimator is

 = N S i=1 1 n n j=1 1 {Yi,j =1} , (11)
where Y i,j ∈ {0, 1} is the number of molecules that are received within a single time slot from the i-th SN in the j-th training slot. Let us define

P0,i = 1 n n j=1 1 {Yi,j =1} , (12)
which is an empirical distribution function. Thus, it follows that the estimator  is consistent [START_REF] Dvoretzky | Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator[END_REF]. Moreover, we observe that

P(|A -Â| > N S ) ≤ P N S i=1 |P 0,i -P0,i | > N S ≤ 1 - N S i=1 1 -P |P 0,i -P0,i | > . (13)
Based on these considerations, we can apply the Dvoretzky-Kiefer-Wolfowitz inequality [START_REF] Dvoretzky | Asymptotic minimax character of the sample distribution function and of the classical multinomial estimator[END_REF], which yields

P(|A -Â| > N S ) ≤ 1 - N S i=1 (1 -2e -2n 2 ) = 1 -(1 -2e -2n 2 ) N S . (14)
Let the time slot duration be T S ∈ {T 1 , . . . , T K }. The bound in [START_REF] Kaelbling | Reinforcement learning: A survey[END_REF] provides a means of obtaining the confidence intervals for the probability of error and the throughput, which is defined as

R l = 1 T l (1 -P e,l ), (15)
for the l-th time slot duration, where P e,l is the probability of error corresponding to a time slot duration T l . In particular, given an estimate Âl it follows that the probability A l ∈ [ Âl -N S , Âl + N S ] is bounded by [START_REF] Kaelbling | Reinforcement learning: A survey[END_REF]. Note that P e,l in ( 9) is a monotonically decreasing function of A, which is due to the fact that the number of molecules that arrive in a single time slot is increasing as A increases.

Using the fact that P e,l is a monotonically decreasing function of A and ( 14), the confidence interval bound on R j

1 -P e,l ( Âl -N S (φ)) T l ≤ R l ≤ 1 -P e,l ( Âl + N S (φ)) T l , (16)
holds with probability greater than 1 -φ, where

(φ) = 1 2n log 2 1 -N S √ 1 -φ . (17)
The confidence interval in (16) bounds the error probability and throughput with n trials, which is useful for providing guarantees on the throughput given a limited training phase. The confidence interval also provides a means for optimizing the time slot duration. Since it is undesirable to obtain tight confidence intervals for every time slot duration, we can use the confidence interval as the basis for an online optimization algorithm. The algorithm can provide guarantees for the probability that the selected time slot duration is suboptimal.

We now present our online algorithm to minimize the probability that a suboptimal time slot duration is selected, as opposed to the time-consuming approach of estimating the throughput for each time slot duration separately. Let 0 < φ < 1 be the confidence level in (16) and the stopping criterion is δ > 0. Suppose that the training phase is in the k-th time slot, then:

1) If 1 ≤ k ≤ K, set j = k and select time slot duration T k .
Using [START_REF] Trees | Detection, Estimation, and Modulation Theory, Part I[END_REF], compute the estimate Âk and the confidence interval via (16). 2) If k > K, select the time slot duration j such that 1 T j

(1 -P e,j ( Âj + N S (φ)))

is maximized. The training phase is ended if for the selected time slot T j , the stopping criterion 1 T j P e,j ( Âj -N S (φ)) -P e,j ( Âj + N S (φ)) < δ (19) is satisfied and the intersection of the confidence intervals defined by ( 16) is empty. Otherwise, the algorithm continues with k ← k + 1.

Note that the algorithm is based on the confidence interval estimate in ( 16), Thus, it follows that the probability that the time slot duration T j selected by the algorithm is suboptimal based on ( 9) is upper bounded by φK.

We remark that the estimation procedure is based on the estimation of the cumulative distribution function for the first passage time. Hence, our algorithm can be applied to systems with any diffusion process as long as the increments are independent and the first passage times for each SN-FN link are each independent.

V. NUMERICAL RESULTS In this section, we show the tradeoffs between performance and complexity of our algorithm by comparing its robustness with time duration spent on the training phase. The system consists of a single FN and N S = 100 transmitting SNs located at equal distance d = 10 -2 m from the FN. To illustrate the behavior of our algorithm in the presence of anomalous diffusion, we consider a fluid medium with a drift velocity v = 0.01 m/s, p E = 0.5, p D = 1, stability exponent α = 1.2 and scale parameter σ = 0.01 for each SN-FN link. We assume that the network can select a time slot duration in the set {1.3, 1.4, 1.5, 1.6, 1.7}.

Fig. 1 compares the throughput for each time slot duration obtained directly from Monte Carlo simulations and from the estimated parameter  using our algorithm for each time slot duration separately with a confidence level φ = 0.01. We observe that the estimate from our algorithm is in good agreement with the throughput obtained from Monte Carlo simulations. This implies that our algorithm can obtain the optimal throughput with probability dependent on the stopping time. Fig. 2 plots the complexity in terms of the number of iterations for the training period as the confidence level varies, which determines the probability the selected throughput is suboptimal. We observe that the stopping criterion has a large influence on the number of iterations. Despite this, we can still obtain an accurate estimate of the throughput in our setting with a small number of iterations as observed in Fig. 1.

VI. CONCLUSION

We have proposed an algorithm to optimize molecular communication in nanomachine sensor networks in the presence of a medium with anomalous diffusion. A key feature of our algorithm is that it does not rely on detailed statistical knowledge of the underlying diffusion processes and that it provides a guarantee on the probability that the throughput is optimal. Simulation results suggest that good performance can be obtained with very short training periods.

Fig. 1 .

 1 Fig. 1. Throughput with different stopping criteria.

3 Fig. 2 .

 32 Fig. 2. Number of iterations with different stopping criteria.