T Le Quang

L Huder

F Lipp Bregolin

A Artaud

H Okuno

S Pouget

N Mollard

G Lapertot

A G M Jansen

F Lefloch

E F C Driessen

C Chapelier

V T Renard

Epitaxial electrical contact to graphene on SiC

Establishing good electrical contacts to nanoscale devices is a major issue for modern technology and contacting 2D materials is no exception to the rule. One-dimensional edge-contacts to graphene were recently shown to outperform surface contacts but the method remains difficult to scale up. We report a resist-free and scalable method to fabricate few graphene layers with electrical contacts in a single growth step. This method derives from the discovery reported here of the growth of few graphene layers on a metallic carbide by thermal annealing of a carbide forming metallic film on SiC in high vacuum. We exploit the combined effect of edge-contact and partially-covalent surface epitaxy between graphene and the metallic carbide to fabricate devices in which low contactresistance and Josephson effect are observed. Implementing this approach could significantly simplify the realization of large-scale graphene circuits.

I. INTRODUCTION

Surface electrical contacts to two-dimensional materials suffer from the poor coupling between the 2D surface and the 3D metal. [1][2][3] The situation is further degraded by contamination in the lithographic processing and/or layer transfer. 4 The best innovation is the recent realization of one-dimensional edge-contacts to graphene 5 possibly combined with large doping. 6 However, the improvement of contact resistance is made at the expense of technological simplicity since the contact fabrication necessitates several steps. Edge bonding [7][8][9] and large electron transfer 10,11 are known to occur during the growth of graphene on SiC and could be exploited for electrical contacts if SiC was replaced by a similar yet conducting material. Conducting carbides appear as good candidates since they have similar chemical properties and since they could allow new functionalities owing to additional material properties such as magnetism or superconductivity. The growth of graphene on carbides other than SiC was first demonstrated by Foster, Long and Strumpf. 12 In 1958, they showed that "aluminum carbide dissociates in the vicinity of 2200-2500 • C, at atmospheric pressures, to aluminum vapor and pure single crystals of graphite" establishing that other carbides could potentially be used for graphene technology. Nevertheless, this subject has remained unexplored owing to the lack of commercial substrates. 13 In this work, we demonstrate that few graphene layer can be grown on a metallic carbide by thermal annealing of a carbide forming metal film (niobium or tantalum) on SiC in high vacuum circum-venting the problem of metallic carbide substrate availability. Based on this discovery we describe a resist-free and scalable method to fabricate few graphene layers (FGL) with electrical contacts in a single growth step. The combined effect of edgecontact 5 and partially-covalent surface epitaxy 14,15 between graphene and the metallic carbide allows us to fabricate devices in which low contact-resistance and Josephson effect are observed.

II. EXPERIMENTAL PROCEDURE

A. Sample preparation

The NbC films were prepared by depositing a 40 nm thick layer of niobium by e-beam evaporation or magnetron sputtering on top of the carbon and silicon terminated surface of 4H-SiC substrates. The samples were then annealed in a RF-induction furnace inside a graphite crucible under a pressure lower than 10 -6 mbar, following a recipe similar to that for graphene growth on SiC 10 . The first annealing step ramps the temperature to 1140 • C in 60 minutes and holds this temperature for 30 minutes for degassing. The next temperature ramp is done in 180 minutes up to 1360 • C and temperature is then kept stable for 18 minutes. The cooling down to room temperature follows a reversed procedure.

B. Characterization

After the thermal treatment, the characteristics of the samples were investigated by several complementary techniques. The morphology was analyzed by scanning electron microscopy (SEM) and transmission electron microscopy (TEM). We used ZEISS Ul-tra+ SEM operated at 5 kV and a probe-aberration corrected FEI THEMIS (S)TEM operated at 200 kV. A protective Ni layer was deposited before specimen preparation for TEM observations by focused ion beam.X-ray diffractograms (XRD) were measured using a PANalytical Empyrean diffractometer equipped with a Cobalt anode (λ Kα1 = 1.7890 Å, λ Kα2 = 1.7929 Å) with an Fe filter to reduce the K β contamination, a Göbel mirror and a 2D-Pixcel detector. For the Raman measurements, a circularly polarized Ar-laser (λ=514.5 nm) at 1 mW of power was used as probe and the scattered light was dispersed by a Jobin-Yvon T64000 spectrometer and collected by a CCD detector. The spatial resolution was better than 1 µm, and the spectral resolution was about 1 cm -1 . Electrical measurements of the superconducting properties of the carbides were performed down to 1 K in a Quantum Design PPMS cryostat using a van der Pauw configuration with a 100 µA drive current. The Josephson junctions measurements were done in a Cryoconcept dilution fridge using phase sensitive detection techniques.

III. RESULTS AND DISCUSSION

A. Formation of the metallic carbide.

We have formed a metallic carbide substrate by vacuum thermal annealing at 1360 • C of a 40 nm thick niobium layer on the carbon face of SiC (See Methods). Figure 1a shows a scanning electron microscope (SEM) image of the Nb film before and after thermal annealing. While the pristine Nb film is flat, the annealed film becomes granular. Since Nb has a very high melting temperature (2477 • C), we expect that this change in morphology is due to a solid state reaction between Nb and the SiC substrate [16][17][18][19] rather than to the melting of the Nb film. The high angle X-ray reflectivity spectra (XRD) shown in Fig. 1b confirms that the film has undergone profound changes. This θ/2θ scan of the annealed sample, reveals that besides the expected peaks of the SiC substrate, one observes peaks at 40.63 • and 87.9 • which we attribute to reflections from crystalline NbC. The relative intensities of the different NbC peaks do not match the expected ones for a randomly oriented polycrystalline phase. They indicate the coexistence of a main [111] textured part with a minor polycrystalline contribution. The φ azimuthal scans of the NbC(311) and SiC(109) (Fig. 1c) reflections show that the textured part of the NbC layer grows with the preferential in-plane orientations NbC(111) [111] parallel to SiC(0001) [1010] and NbC(111)[211] parallel to SiC(0001) [0110]. XRD results are confirmed by TEM analysis which reveal a clean SiC/NbC interface with orientation correlation between NbC and SiC lattices (Fig. 1d). We have observed in most cases NbC [111] oriented parallel to SiC [0001] with several in-plane orientations. XRD indicates a complete transformation of Nb to NbC since neither pure Nb, nor other niobium carbide or silicide were observed. Besides, XRD reveals the presence of cubic silicon carbide which has already been reported at the interface between NbC and SiC. 17

B. Electrical characterization of the carbide

We have measured the van der Pauw resistance of the NbC film studying all possible configurations for current injection and voltage probe contacts 20 . This leads to an accurate estimation of the sheet resistance R s = 3.0 ± 0.3 Ω. One then needs to know the film thickness t to determine the film resistivity ρ = R s t. From the TEM images we find that t = 55 ± 10 nm leading to ρ = 16 ± 5 µΩ.cm. The diffusivity D is obtained independently from the superconducting properties of the film. Figure 2a shows the temperature dependence of the electrical resistivity of the sample for various applied magnetic fields. At zero field, the superconducting transition is sharp and occurs at 12 K. Under external fields, this transition becomes broadened and is shifted to lower temperatures. Figure 2b shows the temperature dependence of the upper critical field defined as the field where the resistivity of the film is about 80% of its value in the normal state. Using the Werthamer, Helfand, Hohenberg dependence of the upper critical field H c2 (T ) for orbital pair breaking 21 (red dashed line in Fig. 2b), we estimate the critical field of our material to H c2 =1.64 T at T=0 K. The diffusivity depends on the value of the slope of H c2 (T ) according to: D = (4k B /πe)(dH c2 /dT ) -1 |T =Tc 22 . We find that the electronic diffusivity of our NbC film is 7.3 cm 2 /s. These transport properties can be related to the good stoichiometry of the film which will be discussed in section III.D. They compete with those of the best NbC reported 23 and the method described here to synthesize NbC could prove useful for the realization of hot electron bolometers where large diffusivity and reasonably large RF impedance is needed. 22 C. Graphene on the metallic carbide.

More interestingly, XRD reveals that the annealing process directly yields FGL at the surface of the NbC layer (XRD peak at 30.5 • ). This is confirmed by the TEM image shown in Fig. 3a where about ten layers of graphene are seen on top of NbC. The small distance between the first graphene layer and NbC (111) surface (0,28 nm, see Fig. 3a) indicates a strong interaction between graphene and NbC. 14 This interaction can be called partially covalent because its strength was shown to be between the weak van der Waals force and the normal covalent bond 14 with strong electron transfer 15 for CVD graphene on metallic carbides. Figure 3b shows that some graphene layers are bonded to a NbC terrace edge.

Figure 3c and d show field TEM views for

different NbC grains covered with FGL details of the connection between NbC and FGL can be seen. The morphology observed Fig. 3d is very similar to that observed on graphene on SiC. 24 On the left of the image (blue box enlarged in Fig. 3e) there are only few graphene layers while on the right of the image (red box enlarged in Fig. 3f) there are much more graphene layers. On SiC, similar morphology was interpreted as resulting from the preferential dissociation of SiC at step edges which leads to nucleation and growth at step edges. This morphology suggests that the growth of graphene on NbC proceeds in a similar way as on SiC 7-9 . However, thermal decomposition of NbC at 1360 • C raises questions since previous studies have revealed that NbC bulk single crystals do not dissociate to graphene by thermal annealing in vacuum at least up to 1800 • C. 14 Another scenario could be that SiC dissociates through the NbC film via diffusion in the bulk or at grain boundaries.

Further experiments demonstrate that the carbidization of Nb and graphene growth occur se-quentially rather than simultaneously. Indeed, no graphene signal is observed on the Raman spectrum of a Nb film annealed at 1140 • C (Fig. 4a). Also, Raman measurements indicate that NbC grown at temperatures below 1360 • C is sub-stoichiometric. Figure 4b shows the low wave number Raman spectrum of three 40 nm Nb films after thermal annealing at 1360 • C, 1140 • C and 1000 • C respectively. These spectra show typical features of NbC 25 . Since NbC crystals inevitably contain C vacancies, quasimomentum conservation is broken and both acoustic (A) and optical (O) phonons over the whole Brillouin zone are therefore Raman active. This results in doublets centered respectively on 170 and 230 cm -1 (A), and 570 and 620 cm -1 (O). We note that we did not find any signature of silicide or other carbide in the Raman response of our samples. Acoustic (A) and optical (O) phonons are seen to shift to lower Raman shift as the annealing temperature increases. This global softening of NbC phonons can be interpreted as a lowering C vacancy concentration 25 . Moreover, additional combinations and overtones around 370 cm -1 (2A,O-A) and 790 cm -1 (O+A) develop when the Nb/C ratio approaches unity 25 . These resonances are observed in the sample an- nealed at 1360 • C confirming our hypothesis that higher temperature annealing allows to saturate the carbide. Comparison with previous results 25 allows to set a lower limit of 0.98 to the carbon/niobium ratio in the sample annealed at 1360 • C. This is further confirmed by the superconducting properties of the NbC films (Fig. 4c). In NbC, T c is extremely sensitive to the carbon content and is degraded from above 12 K for perfect stoichiometry 26 to less than 1.5 K for NbC 0.8 . 27 Figure 4c shows that the resistivity vanishes for a critical temperature of 11.9 K for the sample annealed at 1360 • C while the critical temperature is only 7 K for the sample annealed at 1140 • C. Following Ref. 23 in using the critical temperatures of Ref. 27 as a calibration we get NbC 0.98 and NbC 0.93 respectively, demonstrating that higher annealing temperature allows to saturate the carbide and eventually unlocks the growth of FGL.

E. Fabrication of electrical devices and determination of the contact resistance

Having demonstrated the growth of graphene on a metallic carbide, we have adapted the process to fabricate NbC/FGL/NbC devices where the electrical contact is established during the growth. The overall idea is to structure the initial layer of Nb and exploit the fact that the FGL grow continuously across the junction between SiC and NbC. The main steps of the device fabrication are described in Fig. 5a. Niobium leads are first deposited through a PMMA mask on the substrate at a distance L = 1 -15 µm defining the nominal junction length. The sample is then annealed at 1360 • C in vacuum to form NbC and FGL which covers the entire sample. At the final step of the fabrication, ribbons of variable width W = 3 -10 µm are defined between the electrodes using O 2 plasma. Figure 5b shows an optical microscope image of a SiC chip with 32 devices demonstrating that the technique can be easily scaled up. Figure 5c shows a SEM image of a 5.2×4.4 µm 2 device after the processing and the TEM image of Fig. 5d illustrates that the FGL film is continuous at the interface between NbC and SiC. The FGL film is thicker in the uncovered parts of SiC which further confirms that the growth of graphene on NbC is delayed by the carbidization of niobium.

The devices of Fig. 5b have variable length and width to investigate the NbC/FGL contact resistance by the transfer length method. 28 Figure 5e shows the product of the resistance of the devices and their width R × W as function of their length at room and liquid helium temperatures. Contrary to previously discussed specimens, the devices presented here were grown on the Si-face of the substrate which shows that the growth works similarly on both faces of SiC. As expected, the resistance is independent of temperature and proportional to the length. The extrapolation of the linear dependence to L = 0 µm allows to estimate the specific contact resistance 2R c ≡ R(L = 0) × W . We find R c = 182 Ω.µm which equals the best reported contact resistance to multilayer graphene. 29 We stress here that electrical contacts are realized during the growth which considerably simplifies the sample pro- cessing compared to other methods where the contacts are prepared later on. 29,30 The low contact resistance results from three factors : First, the contact between the metal and graphene is exempt from contamination due to lithography. Indeed, despite resist is used in the initial deposition of Nb, any residues are blown away during the degassing at 1140

• C (See Methods) leaving a clean NbC/graphene interface formed during the growth. Secondly, the peculiar texturation of NbC on SiC leads to a partially covalent epitaxy of graphene on NbC[111]. 14 Large electron transfer specific to this interface 15 was reported which, together with the low graphene/NbC distance lowers the contact resistance. 2 Thirdly, Fig. 3b andd shows that edge contacts 5 are also formed during epitaxy which also improve the contact resistance to multilayer graphene. 29,30 F. Josephson junctions

We investigated the new possible functionalities offered by this technology by realizing Josephson junctions exploiting the superconducting properties of NbC. The junctions were prepared in the same way as described in Figure 5 though here, electron beam lithography was used to achieve shorter channels length. For simplicity, after growth, the FGL were not etched in the form of ribbons between the NbC electrodes. The inset of Fig. 6a shows a SEM image of the junctions for which transport measurements are reported in the main panel. The transport data reported in Figure 6a reveal that a supercurrent flows in this junction below T ≈ 1. of the electrical contacts since Josephson supercurrent is very sensitive to contact quality. The magnetic field dependence of the critical current I c (H) shown in Figure 6b reveals the quantum interference effect between different current paths along the width of the junction. [START_REF] Tinkham | Introduction to Superconductivity[END_REF] The observed dependence I c (H) shows deviations from the usual Fraunhofer pattern due to the peculiar geometry of the junction and subsequent non uniform current distribution. [START_REF] Dynes | [END_REF][START_REF] Barone | Physics and Applications of the Josephson Effect[END_REF][START_REF] Hart | [END_REF][35] Detailed discussions of the transport measurements can be found in supporting discussion 2.

G. Generalization to another metallic carbide

Finally, similar growth and device fabrication can be performed with Tantalum, demonstrating that the method can be generalized to other carbide forming metals. Figure 7a shows a SEM image of a 20 nm thick ebeam evaporated tantalum (Ta) layer af-ter annealing at 1360 • C. The morphology is similar to that of a niobium layer after the same treatment (Fig. 1a) with crystallization suggesting reaction with the substrate. We attribute the dewetting of the layer to the thiner starting layer compared to experiments with NbC. The solid state reaction to form TaC is confirmed by the superconducting temperature of the layer which is 10.15 K, i.e much higher than that of pure tantalum (4.45 K see Ref. 36 ). Similarly to NbC the superconducting transition temperature of TaC is strongly dependent of its stoichiometry and 10.15 K allows to set a lower limit of 0.99 to the carbon/tantalum ratio. 27 The good stoichiometry is confirmed by Raman measurements presented in Fig. 7c where the Raman active optical and acoustic modes and their overtones suggest a low carbon vacancy concentration 25 . Furthermore the Raman measurements indicate that the Ta layer is covered with a graphitic layer as the characteristic G and 2D bands are observed. This set of observations suggests that tantalum follows the exact same route as Nb in reacting with SiC and forming few graphene layers on top.

We verified that the technique previously applied to NbC to fabricate Josephson junctions can be used to fabricate junctions with TaC contacts; Figure 7e shows the current-voltage characteristics of such a junction (SEM image in the inset of Fig. 7) and Fig. 7f the magnetic field dependence. The junction is narrower and the superconducting current therefore set in at higher temperatures. Detailed discussions of the transport measurements can be found in supporting discussion 3.

IV. CONCLUSION

We have demonstrated the formation in a single annealing step of few graphene layers on SiC and their electrical contacts via solid state reactions between a carbide forming metal and the growth of graphene on the metallic carbide and on SiC. This original growth mechanism was exploited to fabricate FGL devices with low electrical contact resistance where the Josephson effect was observed. This contact preparation could be combined with selective graphene ribbon growth on SiC 37,38 to achieve single step circuit fabrication.

Figure 1

 1 Figure 1. a) Scanning electron microscope image of the niobium layer after thermal annealing at 1360 • C. The inset shows a SEM image of the film prior to the annealing. Both scale bars are 200 nm. b) XRD pattern measured on the annealed film. The peaks of NbC are labeled in red, that of graphite in black, those of SiC in blue and those of experimental artifacts in grey. c) φ-scan of the (311) and (109) reflections of respectively NbC and SiC. d) High resolution STEM image of the interface showing that NbC crystal grows in registry with the SiC substrate. The scale bar is 2.5 nm.

Figure 2

 2 Figure 2. a) Temperature dependence of the resistivity for different values of applied perpendicular magnetic field. Bullets represent the temperature for which ρ is 80 % of its normal state value. b) Temperature dependence of the upper critical magnetic field.

Figure 3

 3 Figure 3. a) TEM image of NbC and the graphitic layer that is found on top of NbC. The scale bar is 5 nm. b) TEM image of the edge of a NbC grain showing graphene layers connected to the edge. The scale bar is 5 nm. c) TEM image showing a zoom out of panel a). The scale bar is 20 nm. d) TEM image of another NbC grain covered by FGL. The regions in the blue and red boxes are shown in e) and f). The scale bar is 20 nm in c) and 10 nm in e) and f). Note that a protective Ni layer was deposited before specimen preparation for TEM observations.

Figure 4

 4 Figure 4. a) The black line corresponds to the Raman spectrum measured on a sample annealed at 1140 • C and is compared to that measured on a sample annealed at 1360 • C (red line). While the red spectrum displays typical sp 2 carbon signatures (G, D and 2D), those signatures are absent from the black spectrum indicating that no FGL are formed during the annealing at this temperature. Sharp peaks marked with asterisks are associated to the vibrational modes of O2 and N2 molecules. For clarity, the black spectrum has been shifted vertically. b) Low wavenumber Raman spectra acquired on samples annealed at different temperatures (red line 1360 • C, black line 1140 • C and blue line 1000 • C). c) Low temperature resistivity of NbC films obtained by annealing at 1360 • C (red curve) and at 1140 • C (black curve).

Figure 5

 5 Figure 5. a) The fabrication starts with the deposition of niobium electrodes on bare SiC (left). After high temperature annealing, the metallic electrodes have been transformed to NbC and the entire surface is covered with FGL (center). A ribbon is then etched in the FGL film using O2 plasma (right). b) Optical microscope image of the chip with 32 four-terminal devices. The scale bar is 1 mm. c) Scanning electron microscope image of a 5.2×4.4 µm 2 ribbon. The scale bar is 1 µm. d) Transmission electron microscope image of the junction between NbC and SiC showing the FGL flowing continuously over the two materials. The dashed lines materialize continuous graphene layers at the junction between SiC and NbC. The scale bar is 10 nm. e) Resistance times width of the devices as function of the length. For each length 4 to 5 devices were used to determine the plotted mean value with standard deviation as an error bar. The resistance was measured at room temperature ( ) and at 4 K ( ). The expression R.W = 1400 × L + 364 Ω.µm fits the experimental data (---).

Figure 6 .

 6 Figure 6. Characterization of a Josephson Junction. a) Temperature dependence of the current-voltage characteristic of the junction (T=50 mK, 300 mK, 600 mK, 900 mK, 1.2 K, 1.5 K, 2.1 K ). The inset presents a SEM image of the junction (the scale bar is 200 nm). b) Color plot of the magnetic field and current dependence of the junction's differential resistance at 50 mK.

Figure 7

 7 Figure 7. a) SEM image of a SiC sample covered with an initially 20 nm thick layer of tantalum after annealing at 1360 • C. The scale bar is 400 nm b) Resistance of the TaC film as function of temperature showing a superconducting transition at 10.15 K. c) Low wavenumber Raman spectra acquired on the TaC sample. The shaded area shows the wave numbers where the signal is removed by the spectrometer's filter. d) Raman spectrum obtained on the same sample showing the typical graphene signal (G, 2D and D peaks) e)Temperature dependence of the current-voltage characteristic of a 70 nm junction (T=50 mK, and from 300 mK to 3.3 K in steps of 300 mK). The inset presents a SEM image of the junction (the scale bar is 800 nm). f) Differential resistance of the junction as function of current and magnetic field at 50 mK showing a Fraunhofer pattern.

ACKNOWLEDGEMENT

TLQ was supported by a CIBLE fellowship from Region Rhône-Alpes. FLB acknowledges the support from Zeropova CEA program and Enhanced Eurotalents Marie Curie fellowship. We thank Valérie Reita and the optics and microscopy technological group for valuable support on Raman spec-troscopy, and Nedjma Bendiab for fruitful discussions. We thank C. Marcenat for lending his dilution fridge and technical assistance. We thank F. Gustavo for assistance during metal deposition.

SUPPLEMENTARY DISCUSSION

Supplementary data related to this article can be found below.


Adrien Mazoyer

Rémy Drouilhet

Stéphane Despréaux

Bernard Ycart

flan: An R Package for Inference on Mutation Models

HAL is

Introduction

Mutation models are probabilistic descriptions of the growth of a population of cells, where mutations occur randomly during the process. Data are samples of integers, interpreted as final numbers of mutant cells. These numbers may be coupled with final numbers of cells (mutant and non mutant). The frequent appearance in the data of very large mutant counts, usually called "jackpots", evidences heavy-tailed probability distributions. The parameter of interest is the mutation probability for a mutant cell to appear upon any given cell division, denoted by π. In practice, π is typically of order 10 -9 -10 -11 . Computing robust estimates for π is of crucial importance in medical applications, like cancer tumor relapse or multidrug resistance of Mycobacterium Tuberculosis for instance.

Any mutation model can be interpreted as the result of the three following ingredients:

• a random number of mutations occurring with small probability among a large number of cell divisions. Due to the law of small numbers, the number of mutations approximately follows a Poisson distribution. The expectation of that distribution, denoted by α, is the product of the mutation probability π with the total number of divisions;

• from each mutation, a clone of mutant cells growing for a random time. Due to exponential growth, most mutations occur close to the end of the experiment, and the developing time of a random clone has exponential distribution. The rate of that distribution, denoted by ρ, is the relative fitness, i.e. the ratio of the growth rate of normal cells to that of mutants;

• the number of mutant cells that any clone developing for a given time will produce. The distribution of this number depends on the distribution of division times of mutants.

Using the theory of continuous time branching processes [START_REF] Bellman | On age-dependent binary branching processes[END_REF][START_REF] Athreya | Branching Processes[END_REF], and under specific modeling assumptions, it can be proved that the asymptotic distribution of the final number of mutants has an explicit form. A first mutation model with explicit distribution is the well known Luria-Delbrück model [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF].

Other mathematical models were introduced by Lea and Coulson [START_REF] Lea | The distribution of the number of mutants in bacterial populations[END_REF], followed by Armitage [START_REF] Armitage | The statistical theory of bacterial populations subject to mutation[END_REF] and Bartlett [START_REF] Bartlett | An introduction to stochastic processes, with special reference to methods and applications[END_REF]. In these models, division times of mutant cells were supposed to be exponentially distributed. Thus a clone develops according to a Yule process, and its size at a given time follows a geometric distribution. The distribution of final mutant counts is also explicit when division times are supposed to be constant. This latter model is called Haldane model by Sarkar [START_REF] Sarkar | Haldane's solution of the Luria-Delbrück distribution[END_REF]; an explicit form of the asymptotic distribution is given in Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF]. General division times have been studied by Ycart [START_REF] Ycart | Fluctuation analysis: can estimates be trusted?[END_REF], but no explicit distribution is available apart from the exponential and constant division times.

The first estimation method was given by Luria and Delbrück [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF]. It is based on the simple relation between the probability of null counts in the sample, and the mutation probability, and it is called P0 method. Of course, if the sample does not contain null counts, the method cannot be applied. Apart from the P0 method, all other methods couple the estimation of π or α, with the estimation of ρ. When the distribution of final numbers has an explicit form, the Maximum Likelihood (ML) is an obvious optimal choice [START_REF] Ma | Analysis of the Luria-Delbrück distribution using discrete convolution powers[END_REF][START_REF] Zheng | New algorithms for Luria-Delbrück fluctuation analysis[END_REF]. However, because of the jackpots, likelihood computation can be numerically unstable. There are several ways to reduce tail effects [28, Sec. 2.2], among which "Winsorization" consists in truncating the sample beyond some maximal value. Another estimation method uses the probability generating function (GF) [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF][START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]. The estimators of α and ρ obtained with the GF method proved to be close to optimal efficiency, with a broad range of calculability, a good numerical stability, and a negligible computing time. For the three methods, P0, ML, and GF, the estimators of α and ρ are asymptotically normal. Thus confidence intervals and p-values for hypothesis testing can be computed, for one sample and two sample tests.

The problem with classical mutation models, is that they are based on quite unrealistic assumptions: constant final number of cells [START_REF] Angerer | An explicit representation of the Luria-Delbrück distribution[END_REF][START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF][START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF], no cell deaths (Angerer [2, Sec. 3.1]; Dewanji et al. [START_REF] Dewanji | A generalized Luria-Delbrück model[END_REF], Komarova et al. [START_REF] Komarova | The fixed-size Luria-Delbrück model with a nonzero death rate[END_REF], and Ycart [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF]), fully efficient plating [START_REF] Stewart | Fluctuation analysis: the probability distribution of the number of mutants under different conditions[END_REF][START_REF] Stewart | Fluctuation analysis: the effect of plating efficiency[END_REF][START_REF] Angerer | A note on the evaluation of fluctuation experiments[END_REF], or, as mentioned above, exponential distribution of division times. Using a model for estimation, when the data have been generated by another one, necessarily induces a bias on estimates. For instance, if cell deaths are neglected, mutation probability will be underestimated.

Several informatic tools have already been developed for fluctuation analysis [START_REF] Zheng | Statistical and algorithmic methods for fluctuation analysis with SAL-VADOR as an implementation[END_REF][START_REF] Hall | Fluctuation AnaLysis CalculatOR: a web tool for the determination of mutation rate using Luria-Delbrück fluctuation analysis[END_REF][START_REF] Gillet-Markowska | bz-rates: a web-tool to estimate mutation rates from fluctuation analysis[END_REF]. These tools are quite user-friendly. However, they do not take into account all the possible model assumptions mentionned above. The package flan described here, is dedicated to mutation models, and parameter estimation with the three methods P0, ML, and GF. It includes a set of functions for the distribution of mutant cell counts (dflan, pflan, qflan, rflan) and a graphic function (draw.clone). They treat general models, with fluctuating final numbers, cell deaths, and other division time distributions than exponential and constant. The general estimation function is mutestim. It returns estimates for the parameters α, π and ρ, with the three estimation methods, constant or exponential division times, and cell deaths. As a wrapper, a hypothesis testing function (flan.test) is provided. In order to make the package user-friendly, the functions have been designed to resemble classical R functions, like t.test or rnorm.

The paper is organized as follows. Section 2 is devoted to the probabilistic setting: the hypotheses of the different models are described, and the asymptotic results are explained. In section 3, the three estimation methods are exposed, and the biases described above are discussed. A comparison of the three methods in terms of mean squared errors is provided. The user interface and the Rcpp [START_REF] Eddelbuettel | Seamless R and C++ Integration with Rcpp[END_REF] implementation is treated in section 4; examples of execution are shown in section 5.

Mutation models

In this section, probabilistic mutation models are described. The basic modeling hypotheses are the following:

• at time 0 a homogeneous culture of n0 normal cells is given;

• the lifetime of any normal cell is a random variable with distribution function F ;

• upon completion of the generation time of a normal cell:

with probability π one normal and one mutant cell are produced; with probability 1 -π two normal cells are produced;

• the lifetime of any mutant cell is a random variable with distribution function G;

• upon completion of the lifetime of a mutant cell:

with probability δ the cell dies out; with probability 1 -δ two mutant cells are produced;

• all random variables and events (division times, mutations, and deaths) are mutually independent.

Consider that the initial number n0 tends to infinity, the mutation probability π = πn 0 tends to 0, and the time t = tn 0 at which mutants are counted tends to infinity. The scale of time is supposed to be adjusted so that the exponential growth rate of mutants is 1; thus the exponential growth rate of normal cells is ρ. See Athreya and Ney [4, Chap. IV Sec. 4] or Hamon and Ycart [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF] for the definition of the growth rate (also called "Malthusian parameter"). The expected number of mutations before tn 0 is proportional to n0πn 0 e ρtn 0 , and the asymptotics are assumed to be such that this number converges as n0 tends to infinity to α, positive and finite.

Under the above hypotheses, as n0 tends to +∞, the final number of mutants converges in law to the distribution with PGF:

g(z) = exp (-α(1 -h(z))) , (1)
with

h(z) = ∞ 0 ψ(z, t)ρe -ρt dt , ( 2
)
where ψ(z, t) is the PGF of the number of cells at time t in a mutant clone, starting from a single cell at time 0. Observe that it depends on the lifetime distribution of normal cells F only through ρ. The above result is deduced from the theory of continous time branching processes [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]). The expressions (1) and (2) translate the three ingredients described in the introduction:

1. the Poisson distribution with intensity α models the total number of mutations which occur during the process;

2. the exponential distribution with rate ρ is that of the time during which a random clone develops;

3. the distribution with PGF ψ(•, t) is that of the number of cells in a random clone developing during a time interval of length t. The PGF ψ is the solution of a Bellman-Harris equation [START_REF] Bellman | On age-dependent binary branching processes[END_REF] in terms of δ and G.

Hence the expressions of h as an exponential mixture, and of g as a Poisson compound. In practice, the plating process can be less than 100% efficient. In that case, a random number of mutants will not be counted: if only a proportion ζ of the final population is plated, then each cell will be observed with probability ζ. Denote by Mtot and M the total and the observed numbers of mutants. Given Mtot = m, M follows the binomial distribution with parameters m and ζ. Thus, the PGF g of M is given by:

g(z) = E E z M | Mtot = exp (-α (1 -h(1 -ζ + ζz))) . ( 3
)
The PGF (3) defines a parametrized family of distributions, denoted hereafter by M M (α, ρ, δ, ζ, G) (Mutation Model). This is a family of heavy-tailed distributions, with tail exponent ρ: the higher the fitness, the heavier the tail. This directly influences the number and the amount of jackpots. At this point, the PGF ψ can be given as an explicit expression only for two particular lifetime distributions of mutants: exponential, and Dirac (constant lifetimes). The corresponding mutation models will be denoted respectively by LD(α, ρ, δ, ζ) (Luria-Delbrück) , andH(α, ρ, δ, ζ) (Haldane). The functions dflan, pflan, and qflan compute densities, probabilities, quantiles of LD and H distributions (with ζ = 1).

Assuming that a consistent estimator of α has been defined, the problem in practice is to compute reliable estimates of the mutation probability π. The simplest approach assumes that the final number of cells, denoted by N , is constant. An estimate of π is then obtained by dividing the estimate of α by N . However, even under close experimental monitoring, assuming that the final number of cells is a constant is quite unrealistic. Thus, N must be viewed as a random variable with a certain probability distribution function K on [0, +∞). By analogy with (1), the conditional PGF of the number of mutants given N = n, can be given by the following expression:

g (z | N = n) = exp (-πn(1 -h(z))) .
Or else, the conditional distribution of the number of mutants given N = n is the distribution M M (πn, ρ, δ, ζ, G). The distribution function K is supposed to be known and its Laplace transform is denoted by L:

L(z) = E e -zN = ∞ 0 e -zn dK(n) ,
Thus the PGF of the final number of mutants is given by:

g(z) = ∞ 0 g (z | N = n) dK(n) = L (π(1 -h(z))) . ( 4
)
Remark that if N is constant, (4) reduces to (1) with α = πN . In general, the PGF (4) defines a new parametrized family of mutation distributions, denoted herafter by M M F N (π, ρ, δ, ζ, G, K) (Mutation Models with Fluctuating Numbers of cells).

The two particular cases for the distribution G previously mentioned above (exponential and Dirac) will be denoted by LDF N (α, ρ, δ, ζ, K) (Luria-Delbrück with Fluctuating Number of cells) and HF N (α, ρ, δ, ζ, K) (Haldane with Fluctuating Number of cells). As will be shown in section 3, estimating π by the ratio of an estimate of α by the expectation of N induces a negative bias.

The function rflan outputs samples of pairs (mutant counts-final counts) following M M F N distributions where G is an exponential, Dirac, log-normal or gamma distribution, and K is a lognormal or Dirac distribution.

Statistical inference

Here the three estimation methods P0, ML and GF are described. The main features and the limitations of each method are discussed. The three methods compute estimates of α and ρ, under the LD and H models. When couples (mutant counts-final numbers) are given, estimates of π and ρ are calculated under the LDF N or HF N models.

Even if the probabilities and their derivatives with respect to δ for LD and H distributions can be computed, the variations of the whole distribution as a function of δ are too small to enable estimation in practice (see Ycart [START_REF] Ycart | Fluctuation analysis with cell deaths[END_REF] for more details). Thus, the parameter δ is supposed to be known for the three methods.

In the rest of this section, the three estimators are described, their performances compared in terms of MSE, and the possible sources of biases discussed.

Estimators

P0 estimator:

The first method was introduced by Luria and Delbrück [START_REF] Luria | Mutations of bacteria from virus sensitivity to virus resistance[END_REF] when δ = 0. In that case, the probability of null counts in the sample is e -α . Hence α can be estimated taking the negative logarithm of the relative frequency of zeros among mutant counts. Hence the method cannot be applied if the sample does not contain null counts.

If δ > 0, the probability of null counts in the sample depends also on δ. Assuming δ < 1/2, a fixed point of the PGF ψ(•, t) is the extinction probability of a mutant clone [4, Theorem 1, Chap.I]:

δ * = δ 1 -δ .
By definition, δ * is also a fixed point of the PGF (2). Then the probability of null counts in the sample is e -α(1-δ * ) . A consistent and asymptotically normal estimator of α is given by: α0

= -log (ĝ(δ * )) 1 -δ * , ( 5
)
where ĝ denotes the empirical PGF of the final number of mutants.

Consider now that ζ < 1. Ignoring the the unefficient plating will induce a negative bias. A correction has been proposed by Stewart et al. [26,eq. (41)]. However, it can be used only under model LD(m, 1, 0). Indeed, the general expression of the probability of null counts is e -α(1-h(1-ζ)) , which depends on the fitness ρ. It is still possible to extend the estimator (5) to the case where ζ < 1:

α0 = -log ĝ δ (ζ) * 1 -δ * , (6)
with

δ (ζ) * = δ * -(1 -ζ) ζ .
Remark that ( 6) makes sense only if |δ

(ζ) * | 1.
In particular, if δ = 0, the plating efficiency ζ has to be greater than 0.5. Therefore, the Notice that the P0 method does not directly yield an estimator of ρ. If an estimate is desired, the ML method can be used for ρ only, setting α = α0.

ML estimators:

Since algorithms [START_REF] Embrechts | A limit theorem for tails of discrete infinitely divisible laws with applications to fluctuation theory[END_REF][START_REF] Zheng | New algorithms for Luria-Delbrück fluctuation analysis[END_REF][START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF][START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF] enable to compute the probabilities of the LD, and H models, the ML method seems to be an obvious choice. It can be used on two kinds of samples:

1. sample of mutant counts: In that case, the likelihood is computed with the probabilities of the model LD or H. The parameter of interest is α.

sample of pairs of (mutant counts-final numbers):

In that case, the likelihood is computed with the probabilities of the model LDF N or HF N . The parameter of interest is π.

In both cases, ρ can also be estimated. However, when the sample maximum is large, sums of products of small terms must to be computed [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]. The procedure can be very long and numerically unstable. Thus, the ML estimators can fail for large α and small ρ. In practice, this instability problem is avoided using Winsorization [28, Sec. 2.2], which consists in replacing any value of the sample that pass a certain bound by the bound itself. The bound is 1024 by default, and it could be necessary to increase it. All information above the bound is lost, and in an extreme case where the sample minimum is greater than the bound, irrelevant results will be returned.

In theory, it is also possible to explicit the probabilities of the LD model when ζ < 1. However, these computations have not been done for the H model. Thus the plating process is assumed to be fully efficient when the ML method is used.

GF estimators:

The GF method uses the PGF to estimate the parameter of a compound Poisson distribution [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF][START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF]. Let 0 < z1 < z2 < 1 and z3 in (0 ; 1). The estimators of α and ρ are the following:

αGF (z3) = log (ĝ(z3)) h ρGF (z 1 ,z 2 ) (z3) -1 and ρGF (z1, z2) = f -1 z 1 ,z 2 (ŷ) ,
where ĝ denotes the empirical PGF of the final number of mutants, hx is the PGF (2) with ρ = x, and:

fz 1 ,z 2 (x) = hx(z1) -1 hx(z2) -1 and ŷ = log (ĝ(z1)) log (ĝ(z2)) . ( 7
)
From Rémillard and Theodorescu [START_REF] Rémillard | Inference based on the empirical probability generating function for mixtures of Poisson distributions[END_REF], it can be proved that the couple of estimators ( αGF , ρGF ) is strongly consistent and asymptotically normal, with explicit asymptotic variance [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF].

The GF estimators depend on the three arbitrary values of z1, z2, z3. Those tuning parameters are set to z1 = 0.1, z2 = 0.9, and z3 = 0.8. For more details about the choice of those values, see Hamon and Ycart [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF].

In practice, the GF estimators are quite comparable in precision to ML estimators, with a much broader range of calculability, a better numerical stability, and a negligible computing time, even in the case where the ML method fails. For that reason, we have chosen to initialize the ML optimization by GF estimates, to improve both numerical stability and computing time.

The only practical limitation of this method is the following. A zero of the monotone function fz 1 ,z 2 (ρ) -ŷ must be computed. An upper bound for the domain of research must be given, which can be a problem if the sample does not contain jackpots. However in that case, a mutation model is not adapted.

According to (3), the GF estimators or α and ρ can be extended to the case where ζ 1 as follows:

αGF (z3) = log (ĝ(z3)) h ρGF (z 1 ,z 2 ) (1 -ζ + ζz3) -1 and ρGF (z1, z2) = f -1 1-ζ+ζz 1 ,1-ζ+ζz 2 (ŷ) ,
where fz 1 ,z 2 and ŷ are still given by [START_REF] Boe | Fluctuation analysis of mutations to nalidixic acid resistance in Escherichia Coli[END_REF]. By the same reasoning as Hamon and Ycart [START_REF] Hamon | Statistics for the Luria-Delbrück distribution[END_REF], strong consistency and asymptotic normality of the couple ( αGF , ρGF ) can be proved.

The function mutestim computes estimates and their respective standard deviations for α, π and ρ according to the type of input. Moreover, the estimators mentioned here are asymptotically normal. Thus, one and two sample tests can be performed, using the function flan.test. The null hypothesis will be either fixed theoretical values of α, π, ρ in the one sample case, or a difference of the same in the two sample case.

Comparison of the three estimators

The Figures 1 and2 (drawn using ggplot2) show "maps of usage" of the estimation methods under LD and H models. The methods are compared in terms of the relative MSE of ( α, ρ) defined as:

1 - α α 2 + 1 - ρ ρ 2 . ( 8
)
The RGB code is used: red for GF, green for P0, blue for ML. Twenty values of α between 0.5 and 10, and as many values of ρ from 0.2 to 5, were chosen. Thus 400 couples were considered. For each of them, the following procedure was applied:

1. draw 10 4 samples of size 100 of the LD(α, ρ, 0, 1);

2. for each sample, compute ML, GF and P0 estimates of (α, ρ) under LD model;

3. from the 10 4 estimates, compute the relative MSEs of each method;

4. assign a RGB color according to the MSEs. For each method:

• if the MSE is less than 0.05, assign 1 to the corresponding RGB component;

• if the MSE is greater than 1, assign 0 to the corresponding RGB component;

• else, assign 1 minus the MSE to the corresponding RGB component.

The above experience is also performed considering H model instead of LD. The maps have been drawn with a log 5 -scale for ρ (y-axis). The map can be roughly divided into four distinct parts:

• For (α, ρ) ∈ (0.5 ; 3) × (0.2 ; 2.5), the color is essentially grey: the three methods are more or less equivalent.

• For (α, ρ) ∈ (3 ; 10) × (0.2 ; 3.5), the color is magenta: the ML and GF methods are equivalent.

The P0 method provides estimates with large MSEs or cannot be used because of the absence of null counts.

• For small values of ρ, the color is mainly red: The GF method is the only method with an acceptable MSE. Small values of ρ induce large jackpots. Moreover, the number of jackpots increases with α. Because of the winsorization, the ML and P0 method (which uses ML to estimate ρ) provide estimates with very large MSEs.

• For ρ large, the color is darker and tends to black: the three methods provide estimates with large MSEs, specially for ρ ∈ (3.5 ; 5), where jackpots are very small or absent. In those cases, estimating ρ with the GF method is not possible in practice (see previous sub-section).

Consequently, the GF method will provide a biased estimate for α. The ML method, which [START_REF] Breunig | An almost unbiased estimator of the coefficient of variation[END_REF]. For each of 400 couples of parameters α = 0.5 . . . 10 (x-axis) and ρ = 0.2 . . . 5 (y-axis, log 5 -scale), 10 4 samples of size 100 of the LD(α, ρ, 0, 1) distribution were simulated. The estimates of (α, ρ) were calculated with the three methods. Each method is represented by a color: red for GF, green for P0, blue for ML.

uses the GF estimates to initialize the optimization of the log-likelihood, also provides biased estimates. The P0 method can provide good estimates of α whatever the value of ρ, which explains the presence of green areas at the top of the map. In a case where no jackpots are present in the sample it should be considered that a (heavy tailed) mutation model is not adapted.

Figure 2 is quite similar to Figure 1. There are still two remarkable differences:

• For ρ small and m 2, the three methods seem to be more equivalent under H model than under LD model.

• For ρ large, GF method seems to provide better estimates under H model than under LD model.

The three methods should also be compared in terms of computational time. An illustration on real data will be given in section 5. The slowest method is ML, for the reasons discussed in the previous section. It is even slower when the estimates are calculated under Haldane models H or HF N , when δ is positive, or if the initialization of ρ with the GF method fails. The GF method computes estimates of α and ρ (when possible) in negligible time. The P0 method outputs estimates of α in negligible time, but estimates of ρ are as slow as with ML. [START_REF] Breunig | An almost unbiased estimator of the coefficient of variation[END_REF]. For each of 400 couples of parameters α = 0.5 . . . 10 (x-axis) and ρ = 0.2 . . . 5 (y-axis, log 5 -scale), 10 4 samples of size 100 of the H(α, ρ, 0, 1) distribution were simulated. The estimates of (α, ρ) were calculated with the three methods. Each method is represented by a color: red for GF, green for P0, blue for ML.

Bias evidence

If the model used for the estimation does not correspond to the theoretical model, the estimates can be biased. Four different sources of bias are considered:

1. the final counts are random in the data, constant for the estimation model;

2. cell deaths occur in the data, not in the estimation model;

3. the lifetime distribution is different in the data and the estimation model;

4. the plating process is less than 100% efficient.

In each case, simulation experiments have been made along the following lines:

1. draw 10 4 samples of size 100, under one model;

2. for each sample, compute estimates using another model and the true one if available;

3. compare the empirical distributions of θ/θ, where θ is estimator and θ the true value.

For each figure, red lines mark unit, blue lines mark relative biases of 0.9 and 1.1. According to Figures 1 and2, the GF method is at least equivalent in terms of the relative MSE (8) to the ML and P0 methods. Moreover, it is also the best in terms of computational time. Therefore, the bias evidences will be mainly illustrated with GF method.

Fluctuation of final counts:

When N is constant, the estimate of π is derived by dividing the estimate of α by N . As mentioned in previous section, if N is a random variable, the relation between α and π can be explicit if the distribution K is known. However, this is not the case in practice. Usually, estimates of the expectation and variance of N are available at best. Assume that only the first two moments µ and σ 2 of N are known. Then a first order approximation of the Laplace transform L can be used to reduce the bias. This method is explained in Ycart and Veziris [START_REF] Ycart | Unbiased estimates of mutation rates under fluctuating final counts[END_REF] for the P0 method. It has been adapted to ML and GF estimates. Figure 3 shows the influence of the coefficient of variation C = σ/µ on the ML estimate of π. The estimates were calculated with three different approaches:

• divide ML estimates of α by the empirical mean of N and ignore fluctuations of N (left boxplots);

• directly compute ML with the sample of pairs (mutant counts-final counts) (center boxplots);

• derive from ML estimates of α, taking into account of the empirical fluctuations of N (right boxplots).

According to the visual observations, the bias reduction seems to be working well when either the product πµ or C are small: most of the estimates have a relative bias smaller than 10%. However, the efficiency of the correction decreases as product πµ and C increase. In particular, for larger values of πµ, the bias seems to be smaller without correction. It could be improved with a better approximation of L, that implies knowing or estimating higher moments of N . Another solution is to improve the estimation of C. Here C was estimated by the ratio of the empirical standard deviation of the empirical mean, which is known to be a bad method in terms of MSE [START_REF] Breunig | An almost unbiased estimator of the coefficient of variation[END_REF].

Cell deaths:

The PGFs (1), ( 4) and ( 2) depend on δ. Ignoring cell deaths involves a negative bias on the estimate of α. Assuming the exact value is known, this bias is removed. Figure 4 shows the influence of the death parameter δ on the GF estimate of α. The estimates are calculated with two different approaches:

1. computing GF estimates of α with δ = 0 (left boxplots);

2. computing GF estimates of α with theoretical value of δ (right boxplots).

The visual results show that the negative bias induced by ignoring cells death increases with the value of δ: the relative bias can easily exceed 10% for large values of δ. From the theory of branching processes, the growth process of a mutant clone is supercritical and δ has to be smaller than 0.5. In practice δ is smaller than 0.3. According to the boxplots, the relative bias induced by ignoring cell deaths can reach 0.80. These experiments illustrate also the difficulty to estimate δ. For example, the boxplots at the top right of the figure seems to show that the value of the likelihood for α = 4 and δ = 0 is very close to its value for α = 4 and δ = 0.05.

Lifetime distribution:

As mentioned earlier, the PGF h is explicit only for the LD and H distributions, i.e. when lifetimes are either exponential or constant. This is not the case in practice. If another lifetime distribution is used to simulate the data, and either LD or H are used to estimate the parameters, a bias will be induced on α and ρ. Figure 5 illustrates these observations. It shows the influence of the lifetime distribution on the GF estimates of α and ρ. The samples are drawn assuming the lifetimes are log-normally distributed. The estimates of α and ρ are calculated under LD (left boxplots) and H models (right boxplots).

From the visual observations, the LD and H models can be seen as extreme values for the lifetime distribution:

• both models correctly estimate α;

• the LD model overestimates ρ and has a rather large dispersion of estimated values. The bias seems to increase as α increases;

• the H model correctly estimates ρ. Plating efficiency: Ignoring the plating efficiency induces a negative bias on α and a positive bias on ρ. In practice, the exact value is known. Figure 6 shows the influence of ζ on the GF estimates of α. The estimates are calculated with two different approaches:

1. computing GF estimates with ζ = 1 (left boxplots) ; 2. applying the correction of Stewart et al. [26, eq. (41)] to GF estimates (center boxplots);

computing GF estimates with known value of ζ (right boxplots).

The visual results illustrate first the fact that the correction of Stewart et al. [START_REF] Stewart | Fluctuation analysis: the probability distribution of the number of mutants under different conditions[END_REF] should not be used when ρ = 1: the induced bias on α is negative when ρ > 1, positive when ρ < 1. However, the variance of the estimates obtained with these rectification is smaller than that of the estimates calculated with GF method. Figure 4: GF estimates of α ignoring cell deaths or not. Red horizontal lines mark unit. Blue horizontal lines mark relative bias of 0.9 and 1.1. For each of the 12 sets of parameters α = (0.5, 2, 4, 8) (columns), and δ = (0.05, 0.1, 0.2) (rows), 10 4 samples of size 100 of the LD(α, ρ, δ, 1) distribution were simulated, with ρ = 1. The GF estimates of α and ρ were calculated under LD model. Each boxplot represents the distribution of the 10 4 ratio α/α of the estimates obtained without taking account of cells death (left) and with the theoretical value of δ (right).

Implementation details

The available functions are described here; more details are given in the manual. The behavior of inference functions for inputs which are out of practical limitations is described. Some details about the Rcpp implementation are also provided.

H 0.8 1.0 1.2 1.4 1.6 LD H 0. 1.0 1.2 1.4 LD
H 0.8 1.0 1.2 1.4 1.6 8 1.0 1.2 1.4 LD
Figure 5: GF estimates of α and ρ under LD and H models on data drawn with M M model. Red horizontal lines mark unit. Blue horizontal lines mark relative bias of 0.9 and 1.1. For each of the 12 sets of parameters α = (0.5, 2, 4, 8) (rows), and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the M M (α, ρ, 0, 1, G) distribution were simulated, G being the Log-normal distribution adjusted on Kelly and Rahn's data [START_REF] Kelly | The growth rate of individual bacterial cells[END_REF]. The GF estimates of α and ρ were calculated with the two distributions LD(α, ρ, 0, 1) and H(α, ρ, 0, 1). For each couple (α, ρ), the two first boxplots represent the distribution of the 10 4 ratio α/α obtained by the LD model (left) and the H model (right); the two last boxplots represent the distribution of the 10 4 ratio ρ/ρ obtained by the LD model (left) and the H model (right) user. Those functions have been designed on the principle of the classical distribution functions of R. A graphic function draw.clone is also provided. It represents with a binary tree the growth of a clone starting from a single normal cell with mutation occurrences until a finite time. The function mutestim computes estimates of α or π and ρ, using LD or H models. The three estimation methods are available. Fluctuations of final numbers and cells death are included. It returns estimate(s) of the parameter(s) of interest and the standard deviations. The function flan.test uses asymptotic normality to perform one or two sample hypothesis testing. It has been designed on the principle of the classical hypothesis testing functions of R, such as t.test.

As mentioned in section 3, there are practical limitations for each estimation method. If the inputs of the mutestim function do not respect those limitations, it will output errors or warning messages:

• If δ = 0, the P0 method can not be used if the sample does not contain any null counts. In that case, the mutestim function will return an error message.

• For now, an unefficient plating (i.e. ζ < 1) can be taken into account only with the GF method. If another method is specified, the mutestim function will return a warning message and set ζ at 1. • Issues of the Winzorization parameter w of ML method:

1. If the minimum of the sample is larger than w, then the sample of mutant counts will be constant. 2. If w is too large, then the optimization process can be very long. 3. In the mutestim function w is set at 1024 by default.

• The GF method does not have limitations of usage, even for extreme cases where the ML estimators fail, i.e. samples with theoretical large α and small ρ. However, estimating ρ requires to solve the zero equation discussed in section 3, which is theoretically solvable on R + . In practice the interval of research is bounded. Thus, if the sample does not contain any jackpot, which means that ρ is very large, the zero equation may not have any solution on the interval. In that case, the function will return a warning message, and set the estimate of ρ at 1, and the estimate of the standard deviation at 0. In the mutestim function, the domain of research is [0.01 ; 100].

• Moreover, the initialization of the ML method is done with GF method. Then the domain of optimization is [0.1 × θGF ; 10 × θGF ], where θGF is the GF estimate(s) of the parameter(s) of interest. Then, if the GF method does not succeed to estimate ρ, there is no chance to estimate it with ML. A warning message is returned if the initialization of the estimate of ρ with GF fails.

The function flan.test is a wrapper function of mutestim. It will ouput the same errors or warning messages if its inputs do not respect the practical limitations.

Implementation

Since most functions involve loops that are more expensive in R than in C, flan has been implemented with Rcpp modules. This paradigm provides an easy way to expose C++ functions and classes to R. There are four main classes in the C++ implementation:

• FLAN_Sim: random generation for M M and M M F N distributions. One of its members is a variable of the following type;

• FLAN_SimClone: random generation for clone size distribution according to the lifetimes distribution;

• FLAN_MutationModel: computation of the descriptive functions (probabilities, PGF,...) for LD and H distributions. One of its members is a variable of the following type;

• FLAN_Clone: computation of the descriptive functions for clone size distribution according to the lifetimes distribution.

The Rcpp interface enables also to import into the C++ code any R function. In particular, it is interesting to import the R functions which are already implemented in C. Thus no external C/C++ library is required. The installation remains basic, and the size of the installed package is reduced. For example, the computations of LD distributions involve numeric integrations. The C libraries integration and alglib compute integrals with an accuracy close to machine precision. We could use those libraries but the R function integrate is actually implemented in C. The R function can then be directly called into the C++ code. However, such imports increase the computational time and memory consumption. This is thus unsatisfactory. Another solution uses the package RcppGSL [11, chap. 11], which creates an interface between Rcpp and the C library gsl. Several integration methods are thus available. Since it requires only that gsl is correctly installed, this solution is a good compromise between easy setup and computational cost.

Computing the probabilities for the H distribution with δ > 0 involves squaring high degree polynomials. Such polynomials are easily treated by the package polynom. However its implementation raises memory issues, because of the degree of the polynomials involved. A more efficient way is to use the Fast Fourier Transform. It is provided by the C library fftw3, which can raise some installation issues. The R function fft could be directly called into the C++ code. For the same reasons as for the integrate function, it is more adequate to use the package RcppArmadillo (Eddelbuettel [11,chap.10] ; Eddelbuettel and Sanderson [START_REF] Eddelbuettel | RcppArmadillo: Accelerating R with high-performance C++ linear algebra[END_REF]). This package links Rcpp to the C++ library armadillo, which is dedicated to linear algebra. In particular, it includes a performing Fast Fourier Transform.

Finally, likelihood optimizations in the ML and P0 methods are done with a bounded BFGS optimizer. The package lbfgsb3 provides the eponymous function which is implemented in Fortran. It is much faster than the basic R function optim. werhoff$samples$W11$mc ## [START_REF] Angerer | A note on the evaluation of fluctuation experiments[END_REF] 0 0 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 3 3 4 4 4 4 5 5 This sample does not contain any jackpot, then the theoretical fitness in a mutation model, should be very large. If the GF method is used, it will output a warning message and set the fitness at ρ = 1, as customarily done in the litterature [START_REF] Foster | Methods for Determining Spontaneous Mutation Rates[END_REF]: W <-werhoff$samples$W11$mc # Compute GF estimates of mutations number and fitness mutestim(mc = W, method = "GF") ## Warning in mutestim(mc = W, method = "GF"): Impossible to estimate 'fitness' with 'GF'-method: 'fitness' is set to default value 1 and only mutations number is estimated.

## $mutations ## [1] 0.8792917 ## ## $sd.mutations ## [1] 0.260549 ## ## $fitness ## [1] 1 ## ## $sd.fitness ## [1] 0
Notice that Table 1 of [START_REF] Werngren | Drug susceptible Mycobacterium tuberculosis Beijing genotype does not develop motation-conferred resistance to Rifampin at an elevated rate[END_REF] shows mutants counts and mean final numbers of cells for 1 mL of solution. However, the total volume of the culture is 5 mL. In other words, a plating efficiency of 20% has been applied. The fact is the fitness parameter can not be estimated, even taking into account the plating efficiency: # Volume of each sample: 1 mL # Total volume of each culture: 5 mL # i.e. plating efficiency of 0.2 pef <-0.2 # Compute GF estimates of mutation probability and fitness # taking into account the plating efficiency mutestim(mc = W, plateff = pef, method = "GF") ## Warning in mutestim(mc = W, plateff = pef, method = "GF"): Impossible to estimate 'fitness' with 'GF'-method: 'fitness' is set to default value 1 and only mutations number is estimated.

## $mutations ## [1] 2.804244 ## ## $sd.mutations ## [1] 0.74011 ## ## $fitness ## [1] 1 ## ## $sd.fitness ## [1] 0
Using the P0 method is a way to realize that setting ρ = 1 by default can be misleading. Since this method does not depend on the lifetime distribution, the estimate of α will not depend on the value of ρ.

# Compute P0 estimate of mutations number mutestim(mc = W, fitness = 1, method = "P0")

## $mutations ## [1] 2.525729 ## ## $sd.mutations ## [1] 0.678233
The P0 estimate of α is very different from the GF estimate (when ζ = 1). Consider now the sample of David [9, Tab. 2], which includes rifampin-resistant bacteria counts.

david$D11 ## $mc ## [1] 4 0 1 0 1 0 0 0 0 0 ## ## $fn ## [1] 1.3e+09 9.2e+08 1.3e+09 2.5e+09 1.3e+09 1.6e+09 1.3e+09 2.5e+09 ## [9] 2.5e+09 2.0e+09
Since the 4 value can be seen as a jackpot, the GF method can be used to estimate α and ρ. Now let us compute the ML estimates of π and ρ taking into account or not of the final counts, under the LD model. . The empirical mean of the final counts is denoted by µ, the empirical coefficient of variation by C. Table 1 displays the ML estimates of π and ρ, in the same way as for Figure 3. Their 95% confidence intervals are provided. Comparing the second row to the fourth, one can see that neglecting final number fluctuations induces a bias of order 5% on π, 10% on ρ. From the third row, it turns out that the correction taking into account C, has not improved the estimate of π. Notice also that due to the small size sample, the confidence intervals are quite large. Consider finally the data from Boe et al. [START_REF] Boe | Fluctuation analysis of mutations to nalidixic acid resistance in Escherichia Coli[END_REF]Tab. 4]. The author studied mutations of Escherichia coli from sensitivity to nalidixic acid resistance. 23 samples of resistant bacteria counts are provided. As in the original paper, the 23 samples are concatenated as one. The three estimation methods are compared in terms of computational time on the resulting sample of size 1104. The package microbenchmark is used, evaluating 10 4 times each method on the sample. The methods are compared when only α is estimated (ρ = 1), and when the couple (α,ρ) is estimated. In both cases the estimates are computed under the model LD with δ = 0 and ζ = 1.

Estimates of π

B <-unlist(boeal)

# Concatenation of the 23 samples require(microbenchmark) # Comparing the methods in terms of computational performance # (mutations number only) microbenchmark(P0 = mutestim(mc = B, fitness = 1, method = "P0"), GF = mutestim(mc = B, fitness = 1, method = "GF"), ML = mutestim(mc = B, fitness = 1, method = "ML"), unit = "ms", times = 1e4) The results are shown on Figure 7, as boxplots of timing distributions. Times are in milliseconds and plotted on log-scale.

Estimates of α setting ρ = 1

Estimates of the couple (α,ρ) As mentioned earlier, the ML method is the slowest. The GF method seems to be quite faster

Figure 1 :

 1 Figure 1: Map of usage of the estimation methods under LD model. The map compares the three methods according to their relative MSE[START_REF] Breunig | An almost unbiased estimator of the coefficient of variation[END_REF]. For each of 400 couples of parameters α = 0.5 . . . 10 (x-axis) and ρ = 0.2 . . . 5 (y-axis, log 5 -scale), 10 4 samples of size 100 of the LD(α, ρ, 0, 1) distribution were simulated. The estimates of (α, ρ) were calculated with the three methods. Each method is represented by a color: red for GF, green for P0, blue for ML.

Figure 2 :

 2 Figure 2: Map of usage of the estimation methods under H model. The map compares the three methods according to their relative MSE[START_REF] Breunig | An almost unbiased estimator of the coefficient of variation[END_REF]. For each of 400 couples of parameters α = 0.5 . . . 10 (x-axis) and ρ = 0.2 . . . 5 (y-axis, log 5 -scale), 10 4 samples of size 100 of the H(α, ρ, 0, 1) distribution were simulated. The estimates of (α, ρ) were calculated with the three methods. Each method is represented by a color: red for GF, green for P0, blue for ML.

π = 5 × 10 - 10 π = 2 × 10 - 9 π = 4 × 10 - 9 π = 8 × 10 - 9 CFigure 3 :

 510102109410981093 Figure 3: ML estimates of π ignoring fluctuations of final numbers or not. Red horizontal lines mark unit. Blue horizontal lines mark relative bias of 0.9 and 1.1. For each of the 12 sets of parameters π = (0.5/µ, 2/µ, 4/µ, 8/µ) (columns), and C = (0.2, 0.4, 0.6) (rows), 10 4 samples of size 100 of the LDF N (π, ρ, 0, 1, K) distribution were simulated, with ρ = 1 and K being the Log-normal distribution adjusted to mean µ = 10 9 and coefficient of variation C. The ML estimates of π and ρ were calculated under model LD or LDF N . Each boxplot represents the distribution of the 10 4 ratio π/π obtained with LD model with C = 0 (left), LD model with bias reduction (center), LDF N model (right).

4. 1

 1 User interfaceflan can be split into two distinct parts: the distribution of the final number of mutants, and statistical inference. The functions dflan, pflan, qflan compute densities, probabilities and quantiles of LD and H distributions. The function rflan outputs samples of pairs (mutant counts-final counts) following LDF N , HF N , or M M F N where G is log-normal or gamma distribution and K is lognormal or Dirac distribution. K is adjusted to the mean and coefficient of variation provided by the ρ =

Figure 6 :

 6 Figure 6: GF estimates of α ignoring plating efficiency or not. Red horizontal lines mark unit. Blue horizontal lines mark relative bias of 0.9 and 1.1. For each of the 12 sets of parameters α = (0.5, 2, 4, 8) (rows), and ρ = (0.8, 1, 1.2) (columns), 10 4 samples of size 100 of the LD(α, ρ, 0, ζ) distribution were simulated, with ζ = 0.2 and ζ = 0.05. The GF estimates of α and ρ were calculated under LD model. Each boxplot represents the distribution of the 10 4 ratio α/α obtained with ζ = 1 (left), using the Stewart et al. [26] correction (center), and with the true value of ζ (right).

D

  <-david$D11 Dmfn <-mean(D$fn) # Compute ML estimates and confidence intervals of mutation probbility and fitness # with empirical mean and nul coefficient of variation ft <-flan.test(mc = D$mc, mfn = Dmfn) -10 8.109577 ## attr(,"conf.level") ## [1] 0.95 # Compute ML estimates and confidence intervals of mutation probbility and fitness # with empirical mean and empirical coefficient of variation ft <-flan.test(mc = D$mc, mfn = Dmfn, cvfn = sd(D$fn)/Dmfn) (,"conf.level") ## [1] 0.95 # Compute ML estimates and confidence intervals of mutation probability and fitness # with couples (mc,fn) ft <-flan.test(mc = D$mc, fn = D$fn) (,"conf.level") ## [1] 0.95 The sample of final counts is denoted by D (F N ) 11

Table 1 :

 1 ML estimates of π and ρ of data set from [9, Tab. 2] ignoring fluctuations of final numbers or not. Each row shows the ML estimates of π and ρ and their 95% confidence intervals, deducing from LD model with C = 0 (first row), C = C (second row), and directly with LDF N model (third row).

Figure 7 :

 7 Figure 7: Computational time of the three methods on data from Boe et al. [7, Tab. 4]. Data consist in the 23 samples of boeal concatenated as one. For each method, the estimates of α setting ρ = 1 (left boxplots) or of α and ρ (right boxplots) have been computed under model LD. The timings have been returned with microbenchmark, evaluating 10 4 times each method. Times are in milliseconds and plotted on log-scale.

Examples of usageSome examples on the real data included in flan are provided. Practical limitations, influence of bias sources and comparison of the estimation methods in terms of computational time are illustrated. Consider first the eleventh sample of mutant counts of the werhoff data[START_REF] Werngren | Drug susceptible Mycobacterium tuberculosis Beijing genotype does not develop motation-conferred resistance to Rifampin at an elevated rate[END_REF]:

than the P0 method. However, the estimates of ρ of the P0 method is calculated maximizing the likelihood. This step is more costly in term of computational time. If only α has to be estimated, the P0 method is faster than the GF method.


Zineb Abderrahmane

Gowrishankar Ganesh

Andrea Cherubini

André Crosnier

Zero-Shot Object Recognition Based on Haptic Attributes

Robots operating in household environments need to recognize a variety of objects. Several touch-based object recognition systems have been proposed in the last few years [START_REF] Tanaka | Object manifold learning with action features for active tactile object recognition[END_REF]- [START_REF] Liu | Object recognition using tactile measurements: Kernel sparse coding methods[END_REF]. They map haptic data to object classes using machine learning techniques, and then use the learned mapping to recognize one of the previously encountered objects. The accuracy of these proposed methods depends on the mass of the the training samples available for each object class. On the other hand, haptic data collection is often system (robot) specific and labour intensive. One way to cope with this problem is to use a knowledge transfer based system, that can exploit object relationships to share learned models between objects. However, while knowledge-based systems, such as zero shot learning [START_REF] Palatucci | Zeroshot learning with semantic output codes[END_REF], have been regularly proposed for visual object recognition, a similar system is not available for haptic recognition.

Here we developed [START_REF] Abderrahmane | Zeroshot object recognition based on haptic attributes[END_REF] the first haptic zero-shot learning system that enables a robot to recognize, using haptic exploration alone, objects that it encounters for the first time. Our system first uses the so called Direct Attributes Prediction (DAP) model [START_REF] Lampert | Learning to detect unseen object classes by between-class attribute transfer[END_REF] to train on the semantic representation of objects based on a list of haptic attributes, rather than the object itself. The attributes (including physical properties such as shape, texture, material) constitute an intermediate layer relating objects, and is used for knowledge transfer. Using this layering, our system can predict the attributebased representation of a new (previously non-trained) object and use it to infer its identity.

A. System Overview

An overview of our system is given in Fig. 1. Given distinct training and test data-sets Y and Z, that are described by an attribute basis a, we first associate a binary label a o m to each object o with o ∈ Y ∪ Z and m = 1 . . . M . This results in a binary object-attribute matrix K. For a given attributes list during training, haptic data collected from Y are used to train a binary classifier for each attribute a m . Finally, to classify a test sample x as one of Z objects, x is introduced to each one of the learned attribute classifiers and the output attributes posteriors p(a m | x) are used to predict the corresponding object, provided that the ground truth is available in K.

This extended abstract is a summary of submission [START_REF] Abderrahmane | Zeroshot object recognition based on haptic attributes[END_REF] 1 Interactive Digital Human group IDH, Laboratory for Computer Science , Micro-electronics and Robotics LIRMM, University of Montpellier CNRS, 860 rue Saint Priest, 34090 Montpellier, France.

firstname.lastname@lirmm.fr 2 CNRS-AIST JRL UMI3218/RL 1-1-1 Umezono, 305-8560 Tsukuba, Japan. g.ganesh@aist.go.jp

B. Experimental Setup

To collect haptic data, we use the Shadow anthropomorphic robotic hand equipped with a BioTac multimodal tactile sensor on each fingertip. We developed a force-based grasp controller that enables the hand to enclose an object. The joint encoder readings provides us with information on object shape, while the BioTac sensors provides us with information about objects material, texture and compliance at each fingertip 1 .

In order to find the appropriate list of attributes describing our object set (illustrated in Fig. 2), we used online dictionaries to collect one or multiple textual definitions of each object. From this data, we extracted 11 haptic adjectives, or descriptions that could be "felt" using our robot hand. These adjectives served as our attributes: made of porcelain, made of plastic, made of glass, made of cardboard, made of stainless steel, cylindrical, round, rectangular, concave, has a handle, has a narrow part. We grouped the attributes into material attributes, and shape attributes.

During the training phase, we use the Shadow hand joint readings x sh to train an SVM classifier for each shape, and BioTacs readings x b to train an SVM classifier for each material attribute. SVM training returns a distance s m (x) measure for each sample x that gives how far x lies from the discriminant hyper-plane. We transform this score to an attribute posterior p(a m | x) using a sigmoid function.

C. Zero-Shot Learning

The last step is to classify a test sample x as one of the test objects Z = {z 1 , . . . , z L }. According to the DAP model [START_REF] Lampert | Learning to detect unseen object classes by between-class attribute transfer[END_REF], we introduce x to each attribute classifier to obtain a list of attributes posteriors {p(a 1 | x), . . . , p(a M | x)}. From these, we infer the posterior of a test object z l :

p(z l | x) = a∈{0,1} M p(z l | a)p(a | x) = a∈{0,1} M p(a | z l ) p(z l ) p(a) p(a | x).
(

) 1
Since we have the attributes vector a z l associated to z l , p(a | z l ) = 1 if a = a z l and 0 otherwise. Assuming independence between attributes, we obtain p(a | x) = M m=1 p(a m | x). Replacing in (1) yields:

p(z l | x) = p(z l ) p(a z l ) M m=1 p(a z l m | x). (2)
We use MAP estimator to classify x to the best matching class f (x) = z l :

f (x) = argmax z l p(z l | x). (3)

D. Accounting for Inter-grasp variability in haptic data

Grasps are usually planned for a given object shape. Here we aim for the opposite; to recognize the shape (and material) from a grasp. Our system is thus constrained to work with unplanned rudimentary grasps, which in turn can lead to large variability in the observed haptic attributes. First, it is possible that the size of the object is large and that the robot grasp is not able to enclose the object completely. This can lead to the absence of some attributes in each grasp. Second, the object may be very small, or the grasp orientation improper, such that certain fingers dont touch the object, again leading to variability of haptic information across grasps. Furthermore, even a properly enclosed object may be made of different materials leading to different BioTacs readings from each finger and each grasp.

To cope with these problems, we propose an extension of the DAP model that deals with absent attributes and allows for the integration of information from different fingers. We refer to this extended model as Local DAP, or LDAP. LDAP assumes that if the system reports that an attribute is absent from the object, then the attribute posterior is replaced by a uniform probability to express our uncertainty about the presence of this attribute. Moreover, the posteriors we obtain for a material attribute from the B contacting fingers are averaged to obtain a global posterior of the whole grasp. We can formulate these two modifications in the following formulas, for material attributes:

p(a m = 1 | x) = 1 B B b=1 max(0.5, p(a m = 1 | x b )). (4)
and for shape attributes:

p(a m = 1 | x) = max(0.5, p(a m = 1 | x sh )), (5)
E. Evaluation Fig. 3 gives the confusion matrix when applying LDAP on our test object set. Note again that these objects were distinct from the objects the robot trained on. Our approach shows good performance on most objects. It however does less well with some objects such as the measuring cup and the jar, which was partly expected due to the heterogeneity of the objects itself, and the similarity between the objects.

F. Conclusion

In this work, we designed an attribute-based zero-shot haptic recognition system for recognizing a set of household objects. We used our robotic setup, namely the Shadow hand and BioTac sensors, to collect haptic data from each object, and the DAP model to classify a test sample to one of a set of classes, that the robot has never been trained on.

We are now working to improve and extend the system to better combine information from multiple grasps, consider real-valued attributes [START_REF] Parikh | Relative attributes[END_REF], and to integrate our zero-shot learning architecture with previous state of the art methods for recognizing known objects.

Fig. 1 .

 1 Fig. 1. Direct Attribute Prediction for Zero-Shot Haptic Recognition

Fig. 2 .

 2 Fig. 2. Disjoint training (left) and test (right) object sets.

Fig. 3 .

 3 Fig. 3. Confusion matrix of LDAP

A video of the experiments has been attached to this paper and can be found on the IDH YouTube channel: https://youtu.be/GaGHZXHBDG4