Nicolas Bourgeois
email: nbourgeo@phare.normalesup.orgmarie.cottrell

Marie Cottrell

Benjamin Deruelle
email: benjamin.deruelle@univ-paris1.fr

Stéphane Lamassé
email: stephane.lamasse@univ-paris1.fr

Patrick Letrémy
email: patrick.letremy@univ-paris1.fr

Benjamin Déruelle

Lexical Recount between Factor Analysis and Kohonen Map: Mathematical Vocabulary of Arithmetic in the Vernacular Language of the Late Middle Ages

Keywords: Text Analysis, Factorial Component Analysis, Kohonen Map

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Context

One approach to the understanding of the evolution of science is the study of the evolution of the language used in a given eld. That is why we would like to pay attention to the vernacular texts dealing with practical arithmetic and written for the instruction of merchants: such texts are known since the XIII th century, and from that century onwards and especially after the diusion of the Latin Leonard of Pisa's Liber Abaci, the vernacular language appears more and more as the medium of practical mathematics.

Treaties on arithmetical education were therefore mostly thought and written in local languages. In this process, the XV th century appears as a time of exceptional importance because we can see then how the inheritance of two hundred years of practice transfers into words1 . For the authors of these texts, the purpose was not only to teach merchants but also to develop knowledge in vernacular language, and their books were circulated far beyond the shopkeepers' world, as far as the humanists' circles for example.

1.2

An ob jective of historical research: the study of specialized languages

The work previously done [START_REF] Lamassé | Les traités d'arithmétique médiévale et la constitution d'une langue de spécialité[END_REF]) consisted in the elaboration of a dictionary of the lexical forms found in all the treaties in order to identify the dierent features of the mathematical vernacular language of the time. This done, we have worked on the contexts of some especially important words in order to understand the lexicon in all its complexity, and on the specicities of each text to study the proximities and the dierences between them. In other words, we should like to determine the common language that forms the specialized language beyond the specicities of each text.

2

The data, the ob jectives, the protocol

In order to delimit a coherent corpus among the whole European production of practical calculation education books, we have chosen to pay attention to those treaties which are sometime qualied as commercial (marchand in French) which have been written in French between 1415 and 1520. This last date is the date of the publication of L'arismetique novellement compose of Estienne de La Roche which is closely related to the works of Nicolas Chuquet and which provides in the same time an opening towards the Italian authors, such as Fra Luca Pacioli.

In this way, our corpus is in conformity with the rules of the discourse analysis: homogeneity, contrastiveness and diachronism2 . It contains eight treaties on the same topic, written in the same language and by dierent XV th century authors.

The following table 1 describes some elements of the lexicometric characteristics of the corpus and shows its main quantitative imbalance.

Humanities and Social Sciences traditional protocol

Traditionally on this kind of textual data, HSS researchers use to work on the statistical specicities and on the contextual concordances, since they allow an easy discovery of the major lexical splits within the texts of the corpus while remaining close to the meanings of the dierent forms. clustering methods, combined with co-occurrences analysis -see [START_REF] Martinez | Contribution à une méthodologie de l'analyse des cooccurrences lexicales multiples dans les corpus textuels[END_REF] help us to cluster the texts without breaking the links with semantic analysis. However, such a method of data processing requires a preliminary treatment of the corpus, the lemmatization. It consists in gathering the dierent inected forms of a given word as a single item. It oers the possibility to work at many dierent levels of meaning, depending upon the granularity adopted: forms, lemma, syntax. We can justify this methodological choice here by its eect on the dispersion of the various forms which can be linked to the same lemma, a high degree of dispersion making the comparison between texts more dicult. It must also be remembered that in the case of medieval texts, this dispersion is increased by the lack of orthographic norms. In our case, this process has an important quantitative consequence on the number of forms in the corpus, which declines from 13516 forms to 9463, a reduction of some 30%.

The factorial analysis allowed us to establish a typology of the complete parts of the corpus, based upon all the forms. However, it can be useful to improve this global analysis with a probabilistic calculation for each component of the corpus, by using the table of the under-frequencies [START_REF] Lebart | Statistique textuelle[END_REF]).

It makes it possible to compare the parts of the corpus with each other, taking into account the occurrences of the words and their statistical specicities.

This process has been made with a particular attention to meaning of the word in order to suppress ambiguities : a good example is the French word pouvoir which can be a verb translated by "can" or "may", and which is also a substantive meaning "power".

Finally, to realize a clustering of the manuscripts, we have only kept the 219 words with the highest frequencies. The set of words thus selected for text classication relate both to mathematical aspects, such as operations, numbers and theirs manipulations, as well as to didactic aspects. Their higher frequencies reect the fact that they are the language of the mathematics as they appear to be practiced in these particular texts. Thus, in what follows the data are displayed in a contingency table T with N = 219 rows (the words) and p = 8 columns (the texts) and the entry t i,j is the number of occurrences of word i in text j.

Factorial Correspondence Analysis (FCA)

Factorial Correspondence Analysis is a factorial method which provides the simultaneous representation of both the individuals and their characteristics, that is to say the columns and the rows of a table, in our case the texts (columns) and the words (rows). Figure 1 The rst two factors (43.94% of the total variance) show the diversity of the cultural heritages which have informed the language of these treaties. The rst factor (25.03%) discriminates between the vocabulary according to its relation to the university legacy on the left, and to the tradition of mathematical problems on the right.

On the left, we can observe a group whose strong homogeneity comes from its orientation towards mathematical problems (trouver that is to say "to nd", demander which we can translate as "to ask") and their iteration (item, idem ). That vocabulary can be found most often in both the appendix of Triparty en la science des nombres and Le Traicte de la praticque. Furthermore, there are more verbal forms on this side of the axis than on the other. And we can nd verbs like requerir which means "to require", convenir "to agree", faire "to do", vouloir "to want". Some of them are prescriptive, as devoir "to have to" or vouloir "to want" for example, while others introduce examples, as montrer "to show". All these texts contain a lot of mathematical problems and in a way that texts are more practical.

On the right, the texts of BnF. fr. 1339 and Med. Nantes 456 are clearly more representative of the university culture, containing latin words sequences. They describe basic operations and numbers through words developed around the XII th and the XII th century in the universities, such as digit, article and nombre composé3 .

The second axis (17.91% of the variance) is mostly characterized by the text of BNF. fr. 2050 and also by Kadran aux marchans. It displays words of Italo-Provencal origin, like nombrateur which refers to the division's numerator. Designation of the fraction and operation of division take a signicant part of the information while the most contributory words (for ex. gurer "to draw") allow us to examine another dimension of these works: the graphical representation as a continuation of writing.

Correspondence Analysis displays the particularities of each text, but leaves untouched some more complex elements of the data. For instance, we cannot conclude from this opposition that the authors of the appendix of Triparty en la science des nombres and the Traicte de la praticque are ignorant of the university texts that inspire the other books. Correspondence analysis does not make fully possible the analysis of the attractions. Moreover, we cannot assert that the words which appear in the center of the graph represent a "common vocabulary": as a matter of fact, we ought to analyze all the successive factors in order to build the list of the words constituting the "normal" vocabulary.

SOM algorithm for contingency table

One way to overcome the limitations of the Factorial Correspondence Analysis (FCA) consists of using a variant of the SOM algorithm which deals with the same kind of data, that is a contingency table (see [START_REF] Oja | Kohonen Maps[END_REF] for other applications of SOM to text mining). See [START_REF] Cottrell | Theoretical aspects of the SOM algorithm[END_REF] for a denition of this variant of SOM, we called KORRESP.

The KORRESP algorithm consists in a normalization of the rows and of the columns in order to sum to 1, the denition of an extended data table by associating to each row the most probable column and to each column the most probable row, a simultaneous classication of the rows and of the columns onto a Kohonen map, by using the rows of the extended data table as input for the SOM algorithm.

After convergence of the training step, the modalities of the rows and of the columns are simultaneously classied. In our example, one can see proximities between words, between texts, between words and texts. It is the same goal as in Factorial Correspondence Analysis. The advantage is that it is not necessary to examine several projection planes : the whole information can be read on the Kohonen Map. The drawback is that the algorithm is a stochastic one, and that apparent contradictions between several runs can be troublesome.

In fact, we can use this drawback to improve the interpretation and the analysis of relations between the studied words. Our hypothesis is that the repetitive use of this method can help us to identify words that are strongly attracted/repulsed and ckle pairs.

In its classical presentation [START_REF] Kohonen | Self-Organizing Maps[END_REF], [START_REF] Cottrell | Theoretical aspects of the SOM algorithm[END_REF], the SOM algorithm is an iterative algorithm, which takes as input a dataset x i , i ∈ {1, . . . , N } and computes prototypes m u , u ∈ {1, . . . , U } which dene the map.

We know that self-organization is reached at the end of the algorithm, which implies that close data in the input space have to belong to the same class or to neighboring classes, that is to say that they are projected on the same prototypes or on neighboring prototypes on the map. In what follows we call neighbors data that belong either to the same unit or to two adjacent units. But the reciprocal is not exact : for a given run of the algorithm, two given data can be neighbor on the map, while they are not in the input space. That drawback comes from the fact that there is no perfect t between a two-dimensional map and the data space (except when the intrinsic dimension is exactly 2). Moreover, since the SOM algorithm is a stochastic algorithm, the resulting maps can be dierent from one run to another.

We address the issue of computing a reliability level for the neighboring (or no-neighboring) relations in a SOM map. More precisely, if we consider several runs of the SOM algorithm, for a given size of the map and for a given data set, we observe that most of pairs are almost always neighbor or always not neighbor. But there are also pairs whose associations look random. These pairs are called ckle pairs. This question was addressed by de [START_REF] De Bodt | Statistical tools to assess the reliability of self-organizing maps[END_REF] in a bootstrap frame.

According to their paper, we can dene: N EIGH l i,j = 0 if x i and x j are not neighbor in the l-th run of the algorithm, and N EIGH l i,j = 1 if x i and x j are neighbor in the l-th run of the algorithm, where (x i , x j ) is a given pair of data, l is the number of the observed run of the SOM algorithm.

Then they dene the stability index M i,j as the average of N EIGH i,j over all the runs (l = 1, . . . , L), i. e. M i,j = (1/L) L l=1 N EIGH l i,j . The next step is to compare it to the value it would have if the data x i and x j were neighbor by chance in a completely random way.

So we can use a classical statistical test to check the signicance of the stability index M i,j . Let U be the number of units on the map. If edge eects are not taken into account, the number of units involved in a neighborhood region (as dened here) is 9 in a two-dimensional map. So for a xed pair of data x i and x j , the probability of being neighbor in a random way is equal to 9/U (it is the probability for x j to be a neighbor of x i by chance once the class x i belongs to is determined).

Let Y i,j = L l=1 N EIGH l i,j be the number of times when the data x i and x j are neighbor for L dierent, independent runs. It is easy to see that Y i,j is distributed as a Binomial distribution with parameters L and 9/U . Using the classical approximation of Binomial Distribution by a Gaussian one (L is large and 9/U not too small), we can build the critical region of the test of null hypothesis H 0 "x i and x j are neighbor by chance" against hypothesis H 1 : " the fact that x i and x j are neighbor or not is signicant".

We conclude that the critical region for a test level of 5% based on Y i,j , is

R -[L 9 U -1.96 L 9 U (1 - 9 U ), L 9 U + 1.96 L 9 U (1 - 9 U )]

From ckle pairs to ckle words

We call ckle a word which belongs to a huge number of ckle pairs:

|{i, |M i,j -A| ≤ B}| ≥ T
Unfortunately, it is not quite an easy task to nd an appropriate threshold T . Here we have decided to x it according to data interpretation. The 30 cklest words, whose number of safe neighbors/non-neighbors is between 89 and 119, are displayed in table 2.

contraire "opposite" (89) regle de trois "rule of three" (104) depenser "to expend" (112) doubler "to double" (89) savoir "to know" (105) racine "root" (113) falloir "to need" (93) partie "to divide" (105) chose "thing" (113) meme "same, identical" (93) position "position" (107)

compter "to count" (113) pratique "practical" (94)

exemple "for example" (107) dire "to say" (113) seulement "only" (94) demi "half" (108) nombrer "count" (115) double "double" (97)

garder "to keep"(109) raison "calculation, problem" (116) multiplication (99) science "science" (109) donner "to give" (117) reduire "to reduce" (103)

pouvoir "can" (111) ensemble "together" (117) regle "rule" (103)

se "if" (111) valoir "to be worth" (119) FCA with ckle pairs The combination of both techniques FCA and SOM whose result is displayed in gure 1 is interesting because it preserves properties from the FCA while giving additional information about the center of the projection -which is usually dicult to interpret. Indeed, the identication on the FCA of the ckle forms allows us to control the general interpretation of the factorial graph, where some words nd their place because of the algorithm and not because of their attraction with other forms and with the texts. Remember that, on the rst two factorial axes (see section 2.2), we have observed an opposition between the university legacy, on the right, and a more practical pole with rule, problems and fractions, on the left. It was tempting to support this observation with words such as "practical" or "rule of three". On the other hand, the ckle forms enhancement shows that these words are shared between a lot of dierent texts and not only linked to the treaty of Nicolas Chuquet and the Traicté en la praticque. As a matter of fact, they do belong to all the texts. And we can observe that the rst factor opposes two technical languages that are on either side of a set of common words and these words are obviously not to say necessarily in the center of the FCA (see for instance racine "root").

To conclude, we can see that two levels of interpretation are superimposed: the ckle pairs which reveal the shared lexicon and the factorial map which inserts it in a local interaction system. And because the list is not sensitive to the FCA, we can play on this combination for each successive factorial axis. It is the articulation between these two levels which makes this representation interesting. In the end, the meaning of this new kind of factorial map is quite intuitive and oers easy tools to the argumentation.

Perspectives and conclusion

First, we intend to use the proposed method for other corpus to conrm its capacity to extract specialized vocabulary. In particular we want to make a new experimentation on a corpus of medieval and renaissance prologues of epics. This corpus has been constituted in order to discern the political and ideological appropriations of the chivalric culture in the context of the XV th and XVI th centuries Renaissance. Secondly, the method has to be appropriated by linguists in order to improve it according to their paradigms.

Another challenge will be to work on a statistical characterization of a threshold for the denition of ckle pairs. Indeed, while we managed to dene (through condence intervals) a theoretical frame for reliability of a pair, we still need to infer a similar method for each data.

Finally, we think that we have open a new perspective for clustering through Kohonen maps. Indeed, the study of robust attraction/repulsion between data as well as ckle pairs can be translated into a graph. Then, we can apply methods from graph representation and graph mining in order to get visualization containing more information than a single SOM.

Facteur

  Fig.1. Projection on rst two factors of the FCA. The eight texts appear in frames, the slanted words stand for ckle words (this notion will be dened in section 4)

Facteur

  Fig. 2. Projection on rst two factors of the FCA (zoom on the central part). Only the ckle words are represented.

  Then, the factorial and

	Manuscripts and Title	Date	Author	Number of occurrences	Words	Hapax
	Bibl. nat. Fr. 1339	ca. 1460	anonyme	32077	2335	1229
	Bibl. nat. Fr. 2050	ca. 1460	anonyme	39204	1391	544
	Cesena Bibl. Malest. S -XXVI -6, Traicté de la praticque	1471?	Préhoude? Mathieu	70023	1540	635
	Bibl. des nombres nat. appendix of Triparty en la science Fr. 1346, Commercial	1484	Chuquet Nicolas	60814	2256	948
	Méd. Nantes 456 Bibl. nat. Arsenal 2904,Kadran aux marchans	ca. 1480-90 1485	anonyme Jean Certain	50649 33238	2252 1680	998 714
	Bib. St. Genv. 3143	1471	Jean Adam	16986	1686	895
	Bibl. nat. Fr . Nv. Acq. 10259	ca. 1500	anonyme	25407	1597	730
	Table 1. Corpus of texts and main lexicometric features (Hapax are words appearing
	once in a text).

Table 2 .

 2 30 cklest words among 219 studied

These treaties were written not only in French but also in Italian, Spanish, English and in German.

For further explanations about texts, methodology and purpose of the analysis see[START_REF] Lamassé | Les traités d'arithmétique médiévale et la constitution d'une langue de spécialité[END_REF].

Digit is used for 0 to 9, article for every multiple of ten, and nombre composé is a mixed between article and digit.

Authors Suppressed Due to Excessive Length

Practically, in this study, for each pair of words, we can compute (over 40 experiments) the index M i,j = Y i,j /L, and conclude. Henceforth: if their index is greater than A + B, they are almost always together in a signicant way, the words attract each other. if their index is comprised between A-B and A+B, their proximity depends on the text they belong, they are a ckle pair. if their index is less than A -B, they are almost never neighbor, the words repulse each other.

4 Analysis of ckle pairs and nodes

Identication of ckle pairs

We run KORRESP L times and store the result in a matrix M of size (N + p) × (N + p). The value stored in a given cell i, j is the proportion of maps where i and j are neighbors. White cells stand for pairs that are not neighbors with high probability (proximity happens with frequency less than A -B). Grey cells are not conclusive.


F Cazals

T Dreyfus

A Lhéritier

N Malod-Dognin

PALSE: Python Analysis of Large Scale (Computer) Experiments

Keywords: Experiments, data analysis, statistics, python, scripting PALSE Expériences, analyse de données, statistiques, python, scripting

A tenet of Science is the ability to reproduce the results, and a related issue is the possibility to archive and interpret the raw results of (computer) experiments. This paper presents an elementary python framework addressing this latter goal. Consider a computing pipeline consisting of raw data generation, raw data parsing, and data analysis i.e. graphical and statistical analysis. PALSE addresses these last two steps by leveraging the hierarchical structure of XML documents. More precisely, assume that the raw results of a program are stored in XML format, possibly generated by the serialization mechanism of the boost C++ libraries. For raw data parsing, PALSE imports the raw data as XML documents, and exploits the tree structure of the XML together with the XML Path Language to access and select specic values. For graphical and statistical analysis, PALSE gives direct access to ScienticPython, R, and gnuplot. In a nutshell, PALSE combines standards languages (python, XML, XML Path Language) and tools (Boost serialization, ScienticPython, R, gnuplot) in such a way that once the raw data have been generated, graphical plots and statistical analysis just require a handful of lines of python code. The framework applies to virtually any type of data, and may nd a broad class of applications.

1

On Data Generation and Analysis in Computational Science In computational Science, complex phenomena are modeled and simulated by involved algorithms and programs, which themselves often rely on various parameters. Such a complex investigation strategy raises three core problems. The rst one is the correctness of software, namely its ability to meet specications. This is clearly a non trivial issue, since a number of constraints may be imposed, such as logical constraints, the numerical constraints (numerical correctness is non trivial, since rounding oating point representations according to the IEEE 754 standard may result in erroneous oating point values reported), or performance constraints (requirements on the running time may be mandatory, e.g. for real time applications), etc. The second one relates to the ability to re-generate raw data, a non trivial issue given in particular the numerous parameters that may be involved, and also to the possible random nature of the phenomena investigated. This particular motivated the recent Executable Paper Grand Challenge, proposed by Elsevier. 1 The third one, which is our focus in this note, is concerned with the archival and the analysis of raw data. Handling raw data indeed raises several issues: Raw data perenniality. Make the raw data perennial by ensuring that anyone can make sense out of them even once decoupled from the program that generated them. This is non trivial since the executable and scripts / other programs parsing the output need to co-evolve. Raw data availability. Make the perennial raw data accessible, to allow novel analysis by scientists equipped with dierent methodological tools. Raw data parsing and analysis. Ease the parsing of the raw data, so as to prepare the graphical and statistical analysis.

2

The PALSE Design PALSE is a lightweight python framework addressing the issues just discussed. Its design, which we discuss now, is summarized on Fig. 1.

2.1

Preamble: on Data Formats PALSE assumes that a given (computer) experiments yields one le storing the results. To make these data perennial and foster their availability, PALSE uses XML, since this language is a standard one, and most importantly, provides an abstract structure amenable to high-level querying and ltering operations.

Raw Data Generation

By raw data we refer to the results of some (computer) experimentPALSE can be used to handle data generated by any device, or even manually archived. Archives and their (de-)construction is a tenet of PALSE, and following the Boost library 2 serialization engines. Following the documentation 3 , we use the term "serialization" to mean the reversible deconstruction of an arbitrary set of C++ data structures to a sequence of bytes. Such a system can be used to reconstitute an equivalent structure in another program context. That is, an archive refers to a specic rendering of the aforementioned sequence of byte, and PALSE takes as input XML archives. For computer experiments, the generation of XML data naturally depends on the programming language used, two of them being of particular interest: C++ and python.

In C++, the boost serialization tools accommodate the native C++ types and the data structures of the Standard Template Library Fig. 1). In python, the recursive structure of dictionaries makes conversion of dictionary into an XML tree (and vice-versa) a trivial task.

2.3

Raw Data Parsing and Database Creation Consider a set of XML archives, one per (computer) experiment. PALSE imports each such le as an XML trees. A database corresponds to a set of isomorphic trees corresponding to experiments generated with varying parameters. Therefore, several databases may be used to accommodate several sets of parameters, or to store results from dierent sources. The creation of the trees from the les is delegated to the lxml4 library.

Raw Data Querying

Given the XML trees stored in database(s), the retrieval of the values of interest for graphical and statistical analysis combines features of the XML and XPath languages.

The XML hierarchical representation.

XML provides a hierarchical representation of a document, which may be seen as a rooted ordered tree of nodes called Elements (Fig. 2, and the supplemental Fig. 2). Each Element is characterized by a tag-name, and possibly contains attributes, a text eld and a number of child Elements there is also a optional tail string which is not supported by PALSE):

The tag-name of an Element e corresponds to the name of the start-tag (< tagname >) and the end-tag (< /tagname >) of e in its XML representation. An attribute of an Element e is a variable having a particular value dened in its XML representation with the start-tag of e (< tagname att = ′ value ′ >).

The text eld of an Element e is dened as any text in-between the start-tag and the endtag of e in its XML representation, which is not embedded in-between another (start-tag, end-tag) pair (< tagname > textf ield < /tagname >).

All the children of e are the Elements which are dened in its XML representation inbetween the start-tag and the end-tag of e (< tagname >< childtagname >< /childtagname >< /tagname >).

In the sequel, by data value, we refer either to the value of an attribute or of a text eld. An Element may be also embedded in a namespace, that is represented by a pair (prex, uri). The uri is a string of characters identifying a name or a resource, and the prex is a simple name associated to the uri. Thus, if an Element e having a tag-name tag is embedded in a namespace (prex, uri), its tag-name is represented in the XML le by prex:tag. However, in the XML tree, all prexes are replaced by their corresponding uri, so that the tag name of the Element e is represented in the XML tree by {uri}tag. In the sequel, the term name-value-pair refers to the (tag-name, text-eld) of a specic Element.

The XPath query language. Given an XML document, the XPath language has been developed to address specic parts of this document, see http://www.w3.org/TR/xpath/. We note that an XPath query always specify a path to one or more Element(s). An XPath may be specied in three dierent ways5 : (i) a sequence of tag-names only (e.g tag_1/tag_2), denominating all the Elements having as tag tag_2 and that are children of Elements with tag tag_1.

(ii) or a sequence of tag-names, each possibly enriched by an attribute name (e.g tag_1/tag_2[@att], or tag_1[@att]/tag_2), denominating the subset of the Elements previously described, but having the attribute att.

(iii) or a sequence of tag-names, each possibly enriched by an attribute name and a value of this attribute (e.g tag_1/tag_2[@att = ′ val ′ ], or tag_1[@att = ′ val ′ ]/tag_2), denominating the subset of the Elements previously described, but the attribute att having as value val.

Functions oered in PALSE. The three XPath modes just discussed are used by PALSE functions to query the database, in order to retrieve selected Elements or their data values (either from text elds or attributes). Before describing these functions, we note that each of them operates on the whole database, so that each function returns a list of lists (one per XML tree).

⊲ Function get_all_elements_from_database(xpath_query): targets all the matching Elements. ⊲ Function get_all_data_values_from_database(xpath_query): targets the list of data values of elements: if the XPath query is of type (i), the data value is the value of the text eld; if the XPath query is of type (ii) or (iii), and is ended by an attribute (whose value is given or not), the data value is the value of the attribute att.

To target the rst Element (or data value) only in an XML tree that matches the specied XPath: ⊲ Function get_leftmost_elements_from_database(xpath_query) ⊲ Function get_leftmost_data_values_from_database(xpath_query) Finally, given a list of Elements L, a XPath query X, a data value v and a comparator comp, it is possible to select Elements e of L such that the specied Xpath X from e has a data value v e positively compared with v using comp, a functionality also provided for basic comparators (e.g. on integers):

⊲ Function filter_elements_by_data_values_compare_to(L, X, v, comp) ⊲ Function comparators of integers filter_elements_by_data_values_lower_than_integer(L, X, v)
We also note that upon calling a function returning Elements, the selection can be converted into strings thanks to the function get_data_values_from_elements(L, X), with L and X the arguments described above.

Data Manipulation

The data collected by the previous step typically need to undergo processing before being amenable to analysis. Since the previous step supplies python lists of strings, PALSE provides mechanisms to select, sort, and convert these lists into lists of elementary types. Furthermore, PALSE provides tools for combining lists into dictionaries, for further ltering using regular expressions (for strings) or lists of allowed keys or values.

Data Analysis

The lists just produced are ready for graphical and statistical analysis. PALSE encapsulates functionalities from SciPython and R for computing Pearson, Spearman or Kendall's tau correlations between the collected values, as well as Mann-Whitney rank-sum test. More generally, any package interface with python can be used, e.g. gmpy if unlimited-precision integers / rationals / oats must be used to ensure numerical correctness, gnuplot-py to generate eye-candy plots and charts, etc.

Figure 1 The PALSE workow. The ve steps: (1) the native data types of the application(s) are rst serialized into XML archives, stored into XML les. It is assumed that one le is generated for each (computer) experiment (2) PALSE loads these XML les into databases of XML trees (3) Raw data are extracted from the databases using XPath queries (4) these data are prepared (sorted, ltered, ...) (5) the manipulated data are used to perform various data analyses (statistics, plots, . . . ).

Step 2: raw data parsing

Step 3: raw data querying

Step 4: data manipulation

Step 5: data analysis Database 1

Database 2

F ile1 F ile2 F ile3 F ile4 F ile4 ( , )
Step1 : data generation

Databases of XML trees

Lists and dictionaries

Plots and statistics We illustrate the steps of the PALSE workow (Fig. 1) on the problem of mining the relationship between the number of atoms of a monomeric or large oligomeric protein, and the number of exposed atoms [START_REF] Chothia | Structural invariants in protein folding[END_REF][START_REF] Miller | Interior and surface of monomeric proteins[END_REF], using in-house computer application Intervor [START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF]. The details are provided in the supplemental section 4.1.

(•, •) ( , ) ((•, •, •), (•, •, •)) (♦, ♦, ♦) ((•, •, •), (•, •, •), (•, •, •)) ⊲ Native types, ⊲ Lists, ⊲ Vectors, ⊲ Sets, ⊲ Maps, ⊲ User defined classes . . . (namei, valuei) (namei, valuei) (namei, valuei) (namei, valuei) (namei, valuei) ( , ) (•, •) ( , ) ((•, •, •), (•, •, •)) (♦, ♦) ((•, •, •), (•, •, •))

3.2

Using PALSE to Query Complex Experimental Data Files

As a second example, we exemplify the ability of PALSE to handle data les stemming from an involved experimental process, here les from the Protein Data Bank (http://www.rcsb.org), generated thanks to X ray crystallography or NMR experiments. Such les contain structural information of macro-molecules and macro-molecular complexes, and are notoriously challenging to handle due to the complexity of the information provided (with biological, bio-physical, genetic pieces of information inter-connected).

The code provided in the supplemental section 4.2 shows that PALSE allows dealing with PDB les elegantly. As an example, plotting the distribution of the number of atoms observed in all the (polypeptide) chains of a collection of PDB les merely requires 15 lines of code. The dataset used for this illustration consists of the 82 PDB les, all solved by X ray crystallography, contributed by J. Janin.

4 Supporting Information As a simple illustration, we investigate a simple property of macro-molecular complexes, processing the protein complexes from [START_REF] Lo Conte | The atomic structure of protein-protein recognition sites[END_REF], using the Intervor software [START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF].

More precisely, consider a protein complex involving two partners, say A and B. The surface of a partner consists of the atoms which are exposed (or equivalently not buried) when the partner is alone. These atoms are those which may interact with the other partner. Using Intervor [START_REF] Loriot | Modeling macro-molecular interfaces with intervor[END_REF],

we partition the atoms of a complex as exposed, buried, or interface atoms. The corresponding C++ code and resulting archive are sketched on Figs. 1 and2. Finally, PALSE is used to plot the two statistics of interest (Fig. 4). (Protein complex: PDB code 1aip.pdb). This archive is visualized using google-chrome, which allows inspecting the hierarchical structure. The text circledq in red corresponds to the names of the so-called name-value-pairs; concatenating these strings yields the path used to retrieve values in the tree representation of an XML archive.

4 . 1 Using

 41 PALSE to Handle a Custom Computer Application

Supplemental Figure 1

 1 Palse, step 1: generation of an XML archive. In C++, the creation of the archive consists of calling the boost serialization algorithms for native C++ types and for data structures from the Standard Template Library. Note that each call to the function make_nvp creates a so-called name-value pair. template<class archive> void serialize(archive& ar, const unsigned int version) { using boost::serialization::make_nvp; // serialization for a C++ native type ar & make_nvp("nb_atoms", nb_atoms); ar & make_nvp("nb_buried_atoms", nb_buried_atoms); ar & make_nvp("nb_exposed_atoms", nb_exposed_atoms); ar & make_nvp("nb_interface_atoms", nb_interface_atoms); // serialization for a std::vector, which counts the number // of occurrences of the 20 native amino-acid types ar & make_nvp("AA_count_by_name", residue_count_by_name); } Inria Supplemental Figure 2 Example XML archive generated by the code of Fig. 1.

  The XML trees handled by PALSE A tree T is recursively dened from Elements.

					Element (or node)	Not supported
							by PALSE
		XML tree	start-tag, and optional	text field	children:	end-tag
			attribute with its value			empty iff is a leaf
		= <tag-name attribute-name=attribute-value>	text-field	< /tag-name>	tail-string
	3 Using PALSE
	3.1	Using PALSE to Handle a Custom Computer Application
		Data Structures	Results Files	Databases of XML trees	Lists	Lists and dictionaries	Plots and statistics
		from	C++,
		python, . . .

Inria

Figure

2

See http://www.executablepapers.com/

Boost is a set of peer-reviewed C++ libraries, see http://www.boost.org/

http://www.boost.org/doc/libs/1_49_0/libs/serialization/doc/index.html RR n°8165

lxml is an easy-to-use library for handling XML and HTML in Python, see http://lxml.de/ Inria

http://docs.python.org/2/library/xml.etree.elementtree.html RR n°8165

RR n°8165

Supplemental Figure 3 Palse, step 2 to 5. Note in particular that all the values associated with a specic path in an XML tree are retrieved using the concatenation of these paths, e.g. PDB_le_statistics/nb_atoms to collect the size (in atoms) of all the molecular systems processed.

# Step 2: Create the database of XML trees res_dir = "%s/mols_archives/LoConte_Janin_JMB99/results/palse-example" % getenv_or_die("HOME") Lo_Conte_etal_analysis = PALSE_xml_handle(res_dir) Lo_Conte_etal_analysis.build_Document_Object_Models() This second example illustrates the power of PALSE to handle complex data, to accommodate namespaces, and to perform various ltering steps. As an illustration, we plot an histrogram corresponding to the distribution of the number of atoms found in polypeptide chains of les Inria from the PDB. As seen from Fig. 5, the corresponding python code is minimalist. An example histogram produced from this code is shown on Fig. 6. Supplemental Figure 5 Example: atoms per chains of PDB entries. Step 1. The XML databases of input PDB entries are created using the namespaces required by the PDBML le format (XML format for the PDB). Step 2. Since there are possibly multiple PDB entries, chains and atoms are lists of lists: chains is a list of lists of values of an attribute representing the id of a chain, and atoms is a list of lists of Elements representing an atom_site. Step 3. For each PDB entry, we lter the atoms (we discard heteratoms), and then for each chain of this entry, we count the number of atoms belonging to this chain. Step 4. We plot in a 2D histogram the registered number of atoms per chain using R. Note that the output format is eps. In this section, we sketch the main classes of PALSE, and refer the interested user to the commented source code for more details.

Class Python_dico_vs_XML_Etree. A class performing manipulations between python dictionaries and XML Etrees. This class is meant to (de-)construct XML archives from python dictionaries. All its functions are static members.

Class PALSE_xml_DB. A PALSE database consists of a set isomorphic trees coming from XML archives, which can be queried using the functionality discussed in Sec.


Biran Wang

Pascal Lançon

Céline Bienvenu

Pierre Vierling

Christophe Di

Georges Bossis
email: bossis@unice.fr

A general approach for the microrheology of cancer cells by Atomic Force Microscopy

Keywords: Microrheology, Cancer cell, Atomic force microscopy, Biomechanics and Viscoelasticity

The determination of the viscoelastic properties of cells by atomic force microscopy (AFM) is mainly realized by looking at the relaxation of the force when a constant position of the AFM head is maintained or at the evolution of the indentation when a constant force is maintained. In both cases the analysis rests on the hypothesis that the motion of the probe before the relaxation step is realized in a time which is much smaller than the characteristic relaxation time of the material. In this paper we carry out a more general analysis of the probe motion which contains both the indentation and relaxation steps, allowing a better determination of the rheological parameters. This analysis contains a correction of the Hertz model for large indentation and also the correction due to the finite thickness of the biological material; it can be applied to determine the parameters representing any kind of linear viscoelastic model. This approach is then used to model the rheological behavior of one kind of cancer cell called Hep-G2. For this kind of cell, a power law model does not well describe the low and high frequency modulus contrary to a generalized Maxwell model.

Introduction

The rheological properties of single cells have been investigated in the past using a large variety of methods like micropipette aspiration [START_REF] Sato | Application of the Micropipette Technique to the Measurement of Cultured Porcine Aortic Endothelial Cell Viscoelastic Properties[END_REF], traction or compression between two microplates [START_REF] Thoumine | Time scale dependent viscoelastic and contractile regimes in fibroblasts probed by microplate manipulation[END_REF], optical tweezers [START_REF] Mills | Nonlinear elastic and viscoelastic deformation of the human red blood cell with optical tweezers[END_REF][START_REF] Hénon | A New Determination of the Shear Modulus of the Human Erythrocyte Membrane Using Optical Tweezers[END_REF], magnetic twisting cytometry [START_REF] Wang | Mechanotransduction across the cell surface and through the cytoskeleton[END_REF][START_REF] Maksym | Mechanical properties of cultured human airway smooth muscle cells from 0[END_REF], rotational microrheology [START_REF] Wilhelm | Rotational magnetic endosome microrheology: Viscoelastic architecture inside living cells[END_REF], atomic force microscopy (AFM) [START_REF] Radmacher | Measuring the viscoelastic properties of human platelets with the atomic force microscope[END_REF][START_REF] Lim | Mechanical models for living cells--a review[END_REF][START_REF] Kuznetsova | Atomic force microscopy probing of cell elasticity[END_REF] and many other methods. All these methods (except micropipette aspiration) have in common the tracking of a micro or a nanoparticle either at the surface of the cell or even inside the cell in response to an applied force. Then a model is needed to relate the force-displacement curve to the viscoelastic properties of the cell. If the probe is larger than the characteristic mesh size, , of the actin network (typically =100nm) a continuum approach can be used to model this response. A possible approach rests on a two fluid description with a viscoelastic network coupled by a frictional force to a viscous fluid [START_REF] Levine | Response function of a sphere in a viscoelastic two-fluid medium[END_REF]. The approximation of low Reynolds number can be written as:

  2 Re 1     al G (1)
with a the radius of the probe, l the amplitude of motion of the probe and  the density of the fluid. The order of magnitude of a and l being the micron and the order of magnitude of G() being between 0.1 and 10kPa it is readily seen that this condition holds even for frequencies as high as 100kHz. The other condition needed for the use of a frequency generalized Stokes' law is that a compressive wave in the elastic network would be overdamped on a characteristic length equal to the dimension of the probe; this is realized for frequencies above [START_REF] Levine | Response function of a sphere in a viscoelastic two-fluid medium[END_REF][START_REF] Levine | One-and two-particle microrheology[END_REF]:

  2 2     b a ( 2
)
where  and  are Lamé coefficients and ~/ 2 is the friction coefficient between the fluid phase and the elastic network. Introducing  with G the static shear modulus, the crossover frequency becomes:

2 1              b a (3)
We shall see in the following that the order of magnitude of the relaxation time is the second , so we get  b ~10 -2 ; it means that the deviation from generalized Stokes' Einstein equation is only expected at very low frequencies, that is to say on a long time scale. In AFM experiments, as well as with optical tweezers the bead is not completely immersed in the viscoelastic medium, so this continuum approach must be modified to take into account different boundary conditions. When the motion of the probe is perpendicular to the surface as in AFM or indentation experiments, the basic relation which can be generalized in the frequency or time domain is no longer the Stokes law , but the Hertz law which applies to a spherical probe which settles in an elastic medium Due to the change of the contact surface with time, the time dependent generalization of the Hertz law, which rests on the hypothesis of a linear superposition of delayed response to increments of force or displacement, must be properly written. Such a generalization was described, for instance in Johnson [START_REF] Johnson | Contact mechanics[END_REF] and applied to nanoindentation of cells [START_REF] Srinivasa | An FEM approach into nanoindentation on linear elastic and visco elastic characterization of soft living cells[END_REF]. As the cell thickness is often not very large compared to the indentation depth, a correction of the Hertz formula has been proposed [START_REF] Dimitriadis | Determination of Elastic Moduli of Thin Layers of Soft Material Using the Atomic Force Microscope[END_REF] for elastic media which introduces a function of the indentation depth over the cell thickness.

This correction was used to model the time dependent response to an approximate step displacement, with the approximation that during the relaxation of the force, the indentation depth could be considered as constant [START_REF] Darling | A Thin-Layer Model for Viscoelastic, Stress-Relaxation Testing of Cells Using Atomic Force Microscopy: Do Cell Properties Reflect Metastatic Potential?[END_REF]. Besides the time dependent generalization of the Hertz model including the cell thickness correction, a description of the viscoelastic response function of the material is needed. In principle a deconvolution should allow to extract the creep response function G(t) or the compliance function J(t) but, due the limited range of time -or frequency-in the experiment, it is more direct to fit the experimental curve with a model response function. Furthermore the viscoelastic models can give some insight in the physical process underlying their representation. One of the most common models to describe a viscoelastic solid is the Zener model consisting of a Maxwell element in parallel with a spring representing the zero frequency elastic modulus G 0 , the Maxwell element introduces a relaxation time:

   i i i G (4)
This simple model of viscoelastic solid containing three parameters, can describe quite well the mechanical response of different kinds of cells [START_REF] Darling | Viscoelastic properties of zonal articular chondrocytes measured by atomic force microscopy[END_REF][START_REF] Vadillo-Rodriguez | Dynamic viscoelastic behavior of individual Gram-negative bacterial cells[END_REF], nevertheless some other experiments show a more complex response which can be described by a set of relaxation times or by a power law. In this last case an interpretation of the power law behavior was proposed in terms of a distribution of relaxation times associated with a power law for the distribution of the lengths of the elementary units of the actin network which was supposed to present a fractal structure [START_REF] Balland | Power laws in microrheology experiments on living cells: Comparative analysis and modeling[END_REF]. The mechanical response of a cell to different stimuli is important, because several biological functions are regulated by their contact with the neighboring cells or the extracellular matrix [START_REF] Wang | Control of cytoskeletal mechanics by extracellular matrix, cell shape, and mechanical tension[END_REF]. On the other hand it has recently been demonstrated that a small mechanical stimulus induced by the periodic motion of magnetic nanoparticles at the surface of the membrane can cause the death of the cell [START_REF] Kim | Biofunctionalized magnetic-vortex microdiscs for targeted cancer-cell destruction[END_REF][START_REF] Hu | Nanocomposites with Spatially Separated Functionalities for Combined Imaging and Magnetolytic Therapy[END_REF]. Since a magnetic field can easily be applied in vivo it could be an attractive way to destroy tumors if the magnetic nanoparticles specifically bind to the cancer cells. In this paper our objective is to obtain the rheological characterization of a given type of cancer cells in order to use it in a future work to model the motion of different kinds of magnetic nanoparticles deposited on the surface of the cells. The knowledge of the viscoelastic response of the cell to the AFM probe will allow to determine the deepness at which magnetic particles can penetrate when submitted to an oscillating magnetic field of different frequencies. An atomic force microscope was used and in a first section we describe briefly the operating conditions and the biological material. As explained above, different approximated methods are used to identify the viscoelastic parameters from the experiments, so the second section will be devoted to a presentation of the equation used to deduce the material properties from the indentation depth versus time of a spherical indenter, and the prediction of this equation will be tested against the trajectory of the bead deduced from a finite element simulation which mimics the operating condition of the AFM. In the third section the results are analyzed through different rheological models and it is shown that the generalized Maxwell model gives a better agreement with the experimental data than a power law model.

Materials and methods

Besides imaging, another major application of AFM is force spectroscopy through the direct measurement of tip-sample interaction forces as a function of the indentation depth of the tip into the sample. Compared to conventional rheometers the use of AFM to get rheological properties raises three problems. The first one is the determination of the position of the probe relatively to the surface, the second one is due to the fact that the contact surface between the probe and the biological material is changing during the experiment and at least we have no independent measurement of the position of the probe and of the applied force.

As presented in the introduction there are many papers dealing with these problems and we do not intend to make a review of them here, but just to present what we think to be the best way to overcome them.

Cancer cells

The tumor cell which has been studied is called Hep G2, a human liver carcinoma cell line. HEP-G2 cells were grown at a BD Falcon™ 35mm Easy-Grip™ Cell Culture Dish, in Dulbelcco modified Eagle culture medium (Invitrogen, Carlsbad, CA) containing 10% fetal calf serum, glucose (4.5g/L), glutamine (2mM), penicillin (100units/mL), 10 mg/mL gentamycin in a wet (37°C) and 5% CO 2 /95% air atmosphere. For fluorescence microscopy, HepG2 cells were grown on glass cover slips for 48h and fixed with 4% paraformaldehyde for 15 min. Cells were then permeabilized with PBS 0,1% triton. For microtubule staining, cells were incubated with anti beta-tubulin followed by secondary anti mouse antibody coupled to Alexa-488 (Molecular Probe). For actin staining, cells were incubated with phalloidin (Sigma Aldrich). Nuclei were stained with Hoescht. Cover slips were mounted on slides using Permafluor (Thermo Scientific) and observed with Zeiss Axioimager microscope.

The organization of actin and microtubules are shown separately, respectively on the left and right hand side of Fig. 1. The nucleus is stained in blue. It appears that the microtubule network is homogeneous and dense right above the nucleus.

Force spectroscopy with AFM

Atomic force microscopic experiments were performed with an Agilent series 5100 AFM/STM (Agilent, Santa Clara, CA). The probe of AFM which have been used to obtain the rheology of cancer cells is a spherical borosilicate glass probe (Novascan, Ames, IA). The spring constant of the cantilever is obtained by using an AFM option called Thermal K. which calculates the cantilever spring constant through the use of the equipartition theorem who states that the kinetic energy stored in a system, here on a coordinate which is the deflection of a cantilever from its equilibrium position, is equal to one half of the thermal energy of the system [START_REF] Cook | Practical implementation of dynamic methods for measuring atomic force microscope cantilever spring constants[END_REF]. The spring constant of spherical borosilicate glass probe was 0.047N/m. The radius of the spherical probe was: R= 2.5 m.

The deflection d, of the cantilever is related to the applied force by:

           F t kd k z t h t (5)
Where z(t) is the position of the piezoelectric transducer and h(t) the indentation depth. This equation only applies when the indentation begins that is to say for a given position z 0 (t). For positions above the surface, the deflection can change a little bit for the highest velocity of the cantilever due to the viscosity of the feeding liquid; the corresponding force is subtracted as indicated below. In a "force spectroscopy" mode a first detection of the surface is related to a given indentation force which is taken as a (false) reference surface by the AFM, which needs to be corrected. A fit of the indentation curve with two parameters (the surface position and the shear modulus) is the more usual procedure. Another technique consists in comparing the work done during different indentations which are stopped at a same given trigger force. The work done is the integral of the force displacement curve whose value is quite insensitive to the position of the surface, so, supposing the validity of the Hertz formula, it allows to relate the ratio of the shear modulus at different positions of the cell to the ratio of the indentation depths [START_REF] Ahassan | Relative Microelastic Mapping of Living Cells by Atomic Force Microscopy[END_REF]. This is very convenient for cartography of the shear modulus of the surface of a cell, nevertheless a reference point is needed and, above all, this method relies on the validity of the Hertz equation.

We propose here a different method which allows to determine accurately the contact point of the probe with the surface and the time dependence of the force. The process is exemplified in Fig. 2 where, the upper plot is the applied force versus time, and the lower one is the displacement of the piezoelectric head versus time. The (false) surface is found by setting an arbitrary force for which the piezo will stop and set z=0; this is the initial point at t=0 in Fig. 2 (sketch A), then the piezo retracts during 10s at a given velocity and stops 2µm above this initial position (sketch B) and remains in this position during 10s to reach a plateau (P) on the force curve. Since this plateau in the feeding liquid corresponds to a zero force we use it to find the time t sf and the corresponding position z sf of the piezo where the bead leaves the surface. It only needs to extrapolate the horizontal plateau when the probe is at rest well above the surface (t=15~20s in Fig. 2), to get the time t1 where the bead leaves the surface and then the corresponding position (point 2) z sf of the piezo head which is associated with the surface of the cell: h=0. After that, the cantilever head is lowered at a given velocity (from 100µm/s to 100nm/s), depending on our experiment need, and stopped at a given position(sketch C) which is kept constant during 50s; during this time the probe continues to deep and the force is relaxing(sketch D). Note that the time t 2 where the bead leaves the surface was obtained from the previously determined position of the surface z sf (point 3 left insert) and so the force, F 2 , corresponding to this time t 2 (point 4 on the insert) can be precisely determined Finally, the piezoelectric head is raised in 10s towards its initial position at time t=0; and in all our measurements we recovered the same force as the initial one (red dashed curve), showing that we did not enter into the plastic regime. Finally this procedure also gives us access to the hydrodynamic force on the cantilever: F(t2)-F(t1) which should be subtracted from the total force [START_REF] Alcaraz | Correction of Microrheological Measurements of Soft Samples with Atomic Force Microscopy for the Hydrodynamic Drag on the Cantilever[END_REF]. Even for high velocities, here 3µm in 0.01s, this hydrodynamic force is rather small (around 2nN) compared to the total force and, of course, does not play a role on the relaxation part of the force. The last important thing to note in this Fig. 2 is the fact that when the piezo returns to its initial position (red dash line on both graphs) the initial force is also recovered. It is a proof that there is no plastic deformation of the medium and the recovery time is about 10s in this example.

Microrheology with a spherical indenter

Time dependent Hertz model

The quasi-static force is related to the indentation depth by the Hertz equation [START_REF] Chu | Johnson-Kendall-Roberts Theory Applied to Living Cells[END_REF][START_REF] Guevorkian | Aspiration of Biological Viscoelastic Drops[END_REF]:

32 2 4 3 1    E F R h ( 6
)
where E is the Young modulus and R the radius of the spherical probe. In the absence of information on the Poisson ratio, it is usual to consider the cell as an incompressible medium and to take =0.5 and E=3G where G is the shear modulus; in this case:

3 2 16 3  R F C G h with C (7)
It is interesting to note that the Hertz formula can be recovered, within a constant number, by considering that the average strain, , for an indentation depth, h, is given by:

12 22     h h h a Rh R (8)
with a the radius of the contact area. The average force is then:

1/ 2 3/ 2 22        F G S G Rh G R h (9)
The dependence in R 1/2 h 3/2 predicted by the Hertz law is recovered under the approximation <<1. We shall see in the following that, in our experimental conditions, the Hertz formula needs to be corrected. In normal or cancer cells, more than 98% of molecules are water [START_REF] Freitas | Nanomedicine[END_REF], and their viscoelastic properties result from the interaction between water and a poroelastic medium formed by the actin skeleton of the cell. A quite general description of a viscoelastic solid is given by the generalized Maxwell model consisting in a purely elastic part of shear modulus G 0 in parallel with a series of Maxwell units, representing different relaxation processes modeled by a set of relaxation times [START_REF] Dimitriadis | Determination of Elastic Moduli of Thin Layers of Soft Material Using the Atomic Force Microscope[END_REF][START_REF] Chu | Johnson-Kendall-Roberts Theory Applied to Living Cells[END_REF][START_REF] Guevorkian | Aspiration of Biological Viscoelastic Drops[END_REF].

In the frame of the linear theory of viscoelasticity the increment of stress at a time t resulting from a stepwise increment of strain at time t' is given by:

      ''   d t G t t d t (10)
Where G(t) depends on the viscoelastic model under consideration. For instance, for the generalized Maxwell model, we have:

  12 / // 0 1 2 ... ...           i t tt i G t G G e G e G e ( 11
)
It is important to note that the Eq. 10 is not directly relevant to the case where the contact surface between the medium and the tool is changing with time, which is presently the case. Actually the time dependent generalization of the Hertz model starts with Eq. 7. The increment of force at time t is related to the increment of indentation at time t': dh(t') by:

      32 ''  dF t CG t t dh t (12)
Integrating the Eq.12 gives:

      32 0 ' '' '   t dh t F t C G t t dt dt (13)
or still:

          3 2 3 2 0 0 ' ' '        t F t C G h t G t t h t dt (14)
This is a general formula for a spherical probe indenter, the rheological properties of the material being represented by the stress response function: G(t). The same kind of reasoning with the compliance response function J(t) will give:

          32 0 ' 0 ' ' '        t h t C J F t J t t F t dt ( 15
)
with C'=1/C

The relation between the two response functions is obtained by taking the Laplace transform of the respective convolution products, which will give:

    2  G p J p p (16)
Eqs. 14 and 15 are simply time dependent generalizations of the Hertz formula. If the relative indentation depth h/R, or the indentation depth relatively to the thickness of the sample: h/L, are not much smaller than unity, then these formula need to be corrected, that is what we are going to see in the next section.

Corrections to the Hertz model

The Hertz model is known to be accurate only for small values of the contact area with respect to the radius of the probe. For larger indentation the corrections depend on the Poisson ratio [START_REF] Yoffe | Modified Hertz theory for spherical indentation[END_REF][START_REF] Dintwa | On the accuracy of the Hertz model to describe the normal contact of soft elastic spheres[END_REF]. In order to find a correction adapted to our experimental conditions, we took =0.5 and compared, for a large indentation, the force calculated by FEM with the Hertz equation for a purely elastic cylindrical plate whose thickness and radius were very large compared to the probe radius (R=2.5m). The thickness L and the radius of the cylinder were both taken equal to 200µm. The Young modulus E was 575Pa corresponding to a value measured on a cancer cell. We found a small difference (4% for h=1m) between the Hertz equation and the result obtained by FEM. Introducing a first order correction factor in the Hertz equation (Eq.17), gave a result which very well reproduces the FEM result until h/R=0.4 :

  3/ 2 2 4 1 10
31       E R h Fh R (17)
On the other hand the Hertz equation was corrected for the effect of the finite thickness, L, of the sample by [START_REF] Srinivasa | An FEM approach into nanoindentation on linear elastic and visco elastic characterization of soft living cells[END_REF]:

    3/ 2 2 4 31    ER F h f h ( 18
)
where:

  2 2 3/ 2 2 2 3 3/ 2 3 2 0 0 0 0 0 0 0 2 2 3 3 4 4 2 4 16 8 4 3 1 15 5                                   R R R R f h h h h h L L L L (19) 0 2 1.2876 1.4678 1.3442 1     vv v 0 2 0.6387 1.0277 1.5164 1     vv v
Once again, there is a small difference between this modified Hertz equation Eq.18 and FEM result, but this difference disappears if we use the same correction as in Eq.17 Finally the corrected expression of the Hertz model, Eq.20, is found to be in excellent agreement with the FEM result at least for indentation below 20% of the plate thickness.

    3/ 2 2 4 1 10 31       E R h F h f h R (20)
The time dependent generalization Eq.20 is straightforward. Let us call g(h) the function which multiplies the Hertz equation:

    1 10     h g h f h R (21)
Then, instead of applying the time superposition to h 3/2 it should be applied to h 3/2 g(h) so Eq.13 becomes:

          3/ 2 0 ' ' ' ' '     t d F t C G t t h t g h t dt dt (22) or                   3/ 2 3/ 2 0 0 ' ' ' '          t F t C G h t g h t G t t h t g h t dt ( 23
)
Before applying Eq. 23 to find the function G(t) from an experiment it is worth testing its application to a viscoelastic solid represented by the popular Zener model.

Derivation of the probe trajectory in a viscoelastic medium

We wish to predict the trajectory of the spherical probe inside a viscoelastic solid represented by the Zener model consisting of a spring of modulus G 0 , in parallel with a Maxwell element represented by a spring G 1 and a dashpot of viscosity  1 . For this model the response function to a unit stress step is:

  1 01    t z G t G G e (24)
where 

=

1 /G 1 is the relaxation time of the Maxwell branch. Then from Eq.23

            3/ 2 0 1 1      F t C G G h t g h t F t (25) with         1 ' 3/ 2 1 1 1 0 ' ' '       tt t

G F t e h t g h t dt

Deriving Eq. 25 and taking into account that:

          3/ 2 1 1 1 11     Ft G F t h t g h t (26)
and (cf. Eq.5) that:      

   F t k z t h t

We shall end up in the case of the Zener model with the following differential equation:

              3/ 2 0 1 1 1 3/
                               τ τ τ dz h z h G hh C h f h dt R dh dt h h h G G Q h f h h f h R R R (27)
where  

                               R R R R Q h h h h h L L L L ( 28
)
In order to check the validity of this approach the indentation and the force versus time are calculated from Eq. 27 and compared with the result obtained from a simulation done with the software Abaqus which contains the possibility to introduce the viscoelastic parameters of the Zener model. A 2d axisymmetric geometry was used; the cantilever was represented by a spring of stiffness 0.05N/m and the spherical probe by a rigid sphere (Fig. 3). The cell was modeled by a cylinder of thickness L whose viscoelastic properties, defined by the Zener model, were the following: G 0 =191.6Pa, G 1 =118.4Pa,  1 =8.4s which are values obtained for a given cancer cell.

Fig. 3: Calculation of the force with respect to the indentation in a viscoelastic medium with a spherical probe mounted on a spring using the software Abaqus

The result is shown in Fig. 4 where the force was plotted versus time when a displacement at a velocity of 0.1m/s during 10 seconds was applied to the top of the spring and then stopped at z=1µm. The numerical solution of the differential equation Eq.27 was compared to the Abaqus result for a cylinder whose thickness is infinite (in practice L=200µm in Abaqus) and for finite thicknesses of 5 and 10µm. Fig. 4 shows a very good agreement between FEM calculation and the solution of Eq.27 in all the cases, so we can trust the validity of the corrections applied to the Hertz formula and its generalization to time dependent forces. It is also worth noting that it applies for any function z(t) and in particular for a sinusoidal one. Another interesting point revealed by the FEM calculation is that the maximum stress generated by the motion of the spherical probe is not situated directly on the surface of the probe but below, on the axis of indentation (cf. Fig. 3) This maximum stress depends on the indentation depth and, as long as no plastic flow is experimentally observed, it determines the range of stress that the cell can support without yielding. In Fig. 5 is plotted, in the case of a purely elastic medium, the evolution of the maximum stress with the indentation depth. The average stress can be approximated by putting: =F/(2Rh) with F given by Eq.20. Of course this average stress is lower than the maximum stress (cf. Fig. 5), but it is quite remarkable that the ratio between the maximum stress and the average one, remains approximately constant with the indentation at a value of about 1.6 up to 0.6m and then increases slowly until 1.8 for a deepness of 1.7µm.

Deriving the material properties

Eq.25 is a general equation relating the force, F(t), and the indentation, h(t), which are both known experimentally, to the response function G(t). In principle it is possible to extract G(t) by taking the Laplace transform of the convolution product in order to get G(p) and then to come back by inverse Laplace transform to G(t). Nevertheless ,the time scale on which the force is measured is not infinite as in a convolution product and the experimental data can be noisy, so the error coming from this double Laplace transform is difficult to evaluate; furthermore , even if we can get a reliable result for G(t), we would like to compare it with known models. In practice it is much easier to directly begin with a given model for G(t)-like for instance the generalized Maxwell model-and to see afterwards if this model can well represent the experimental curve. As an example of this procedure let us suppose that the function G(t) of Eq. 11 well describes the viscoelastic properties of the material; then inserting it in Eq.23 gives:

        3/ 2 01 ()       nn ii ii F t C G h t g h t F t ( 29
)
where

      1 ' 3/ 2 0 ' ' '         tt t i i i

G F t e h t g h t dt

This is just the generalization of Eq.25 to a set of relaxation times.

In a fit of the parameters G 0 , G i ,  i , (i=1, n), the force F(t) is the experimental function whose values are known for a given set of time values: F exp (t j ) (j=1,N) and the left-hand side of Eq.( 29) is the "theoretical" part: F th (t j ,G 0 ,G i , i ) containing the unknown parameters. The minimization of the residue, using the Levenberg-Marquardt method allows to find the unknown parameters. This procedure can, of course, be used for any kind of function G(t) and will be used in the following both for the generalized Maxwell model (Eq.11) and a power law relaxation function: G(t)=At 

Results

As explained in section 2.2 all the experimental curves have been obtained in force spectroscopy mode with different velocities of the piezoelectric transducer on which is attached the cantilever head. The final position was corresponding to an indentation depth of about 1m. The small viscous force coming from the hydrodynamic friction of the cantilever in the feeding liquid was subtracted and the indentation depth, h(t), was obtained as shown in Fig. 2. The indentation was stopped at a time t=t s . If there was no relaxation during the rising time, or in other words if t s << i , then the terms F i (t) in Eq.29 are negligible and the force at t=t s can be expressed as:

        3/ 2 0       n s i s s i F t C G h t g h t ( 30
)
This is the hypothesis which is usually done [START_REF] Darling | Viscoelastic properties of zonal articular chondrocytes measured by atomic force microscopy[END_REF][18] when dealing with the relaxation part for t>>t s . It means that at t=t s , which is the starting point of the relaxation curve, Eq.30 gives access to the sum of the shear moduli of the model, which is called the high frequency modulus:

0 1 2      G G G G ( 31
)
On the contrary, if the time t s is longer than all the relaxation times, the integrand representing F i (t) in Eq.29 will be different from zero only for t-t'< where t s , so the term containing h(t) can be taken out of the integral in F i (t) with its final value h(t s ); then considering that t s >> i we shall end up with:

        3/ 2 0  s s s F t CG h t g h t ( 32
)
This is now the low frequency limit of the shear modulus which is measured at t=t s but in this situation the indentation is a quasi static one and after stopping there will be no relaxation. This is the way to obtain the static shear modulus G 0 . All the experiments that we are going to discuss were made on the central part of the cell where the height of the cell is maximum and is about 4 m. Actually, as shown in Fig. and in the table just below, the values of G 0 do not vary significantly for different points in the region above the nuclei in the white zone (points E, F, D), but they increase when the measurement is done closer to the border of the cell. These measurements were made in a quasi static way; on the other hand, if only the relaxation part of the force is used then t s needs to be as small as possible to minimize the error, but in any event it will be safer to fit the full curve from the beginning of the indentation. In order to illustrate this point we have presented in Fig. 7 the values of the shear modulus which are obtained from Eqs. [START_REF] Dintwa | On the accuracy of the Hertz model to describe the normal contact of soft elastic spheres[END_REF] or [START_REF] Wilhelm | Elasticity of Stiff Polymer Networks[END_REF] after an indentation made at different velocities on the same cell. It clearly shows that the shear modulus obtained in this way depends on the indentation velocity but this value has no clear meaning in the frame of the viscoelastic models and only a fit of the total indentation curve can allow to attribute safely the values of the fitted parameters to a given model as we shall see in the next section.

Response function G(t) of the cells Hep G2

A typical curve representing the force versus time for a "slow" velocity of the cantilever: v=0.25µm/s is plotted in Fig. 8; the indentation depth corresponding to the maximum of the force is h s =1.08µm. The force (black square line) passes through a maximum when the cantilever stops and then decreases as the probe continue to go down until the elastic stress becomes equal to the stress applied by the cantilever. This curve was fitted with the generalized Maxwell model with a single branch (Zener model) for the 3 parameters: G 0 , G 1 ,  1 .

Fig. 8: Applied force versus time during an indentation with a velocity: v=0.25µm/s and a maximum indentation depth: h s =1.08µm .The blue triangles represent a fit with a power law for G(t) and the red circles a fit with a one branch Maxwell model

As can be seen in Fig. 8, the fit is good on all the time scale for the following parameters:

G 0 =198Pa, G 1 =101Pa,  1 =7
.1s. This model was also observed to give a good representation of the viscoelasticity of some bacterial cells [START_REF] Vadillo-Rodriguez | Dynamic viscoelastic behavior of individual Gram-negative bacterial cells[END_REF] and of chondrocytes [START_REF] Darling | Viscoelastic properties of zonal articular chondrocytes measured by atomic force microscopy[END_REF]. Nevertheless it is worth noting that in these previous works the fit was only done on the relaxation part of the curve, supposing that the indentation was fast enough so that the medium did not have time to relax during the indentation. If this hypothesis was used for the experiment represented in Fig. 8 the relaxation part would be fitted by the function:

    1 3/ 2 01 when             s tt s s s s F t t Ch g h G G e t t ( 33
)
where t s and h s are respectively the time and the indentation corresponding to the maximum of the force. In our example it would give G 0 =203Pa, G 1 =70Pa,  1 =8.4s instead of G 0 =198Pa, G 1 =101Pa,  1 =7.14s. The error is not that big, the maximum error being 40% for G 1 , nevertheless the degree of approximation of Eq.33 depends on the importance of G 1 compared to G 0 and on the relaxation time  1 compared to t s since it supposes that for t=t s the modulus is equal to G 0 +G 1 , that is to say that the material did not have time to relax during the rising indentation.. On the other hand it is clear that the velocity of the cantilever head will give a limit on the smaller relaxation time which can be detected: roughly speaking a relaxation time which would be smaller than t s would not be detected. In practice our Hep G 2 cell gives a good example of this remark. In Fig. 9 is plotted the force versus indentation for a velocity of the piezoelectric head: v=250µm/s and about the same indentation depth h s =1.27µm. The maximum of the force is higher than in Fig. 8, and above all, shows a very fast relaxation of large amplitude. A fit of this curve with a single relaxation time,  1 , (red circle curve in Fig. 9) no longer works, and it is necessary to add another relaxation time,  2 , to give a satisfying agreement (blue triangle curve). The results obtained on several experiments both for fast and slow indentation velocities are summarized in Tab 1. In parentheses are the values calculated without taking into account the thickness of the cell which overestimates the moduli of about 15-20% because the increase of stiffness due to the presence of the bottom wall was attributed to the elasticity of the cell. The rapid relaxation time  2 =0.06s associated with a modulus G 2 of about 450Pa means that some elastic part of the cytoplasm, whose elasticity is more than two times larger than the static one, is associated with a rather low friction during its relaxation ( 2 = 2 G 2 =27Pa.s) compared to the one associated with the slow relaxation process ( 1 = 1 G 1 =450Pa.s). In previous experiments, either with AFM [START_REF] Vadillo-Rodriguez | Dynamic viscoelastic behavior of individual Gram-negative bacterial cells[END_REF] or with specific devices like microplates [START_REF] Thoumine | Time scale dependent viscoelastic and contractile regimes in fibroblasts probed by microplate manipulation[END_REF], the time scale of the indentation or of compression in the case of microplates was much longer so a Zener model with a single Maxwell time as in Fig. 8 was able to represent the viscoelastic behavior of the cell and this fast relaxation time associated with a lower viscosity was not detected. The first Maxwell branch could represent the entropic part associated with the shear elongation of the crosslinks and the second Maxwell branch with a faster relaxation and a higher modulus could come from the change of entropy associated with the compression of the main filaments of the network [START_REF] Wilhelm | Elasticity of Stiff Polymer Networks[END_REF]. Of course many others explanations can be proposed, but the interest of the Maxwell model is that it can be interpreted more easily than other models on a physical basis. Another model which is often used to represent the viscoelastic behavior of gels consists in a power law distribution of the relaxation times; it was also argued to represent the viscoelastic response of the actin network, leading to a power law for the creep relaxation function [START_REF] Balland | Power laws in microrheology experiments on living cells: Comparative analysis and modeling[END_REF]:

Slow indentation G

    J t At ( 34
)
The Laplace transform of J(t) is:

    1 1      J p A p (35)
Since the Laplace transforms of G(t) and J(t) are related by G(p)J(p)=p 2 (cf. Eq.16) it is easy to show that G(t) will be given by:

          0 0 tt 11 , 0 1 11 tt                B Gt A (36) Where     1 11       B A (37)
Formally we should have t 0 =0 but this time shift is an ersatz which is used to avoid the divergence of G(t) in t=0 [START_REF] Schiffmann | Nanoindentation creep and stress relaxation tests of polycarbonate[END_REF]. In practice the fit is realized on the three parameters, B,,t 0 , knowing that t 0 should be small compared to the relevant relaxation times. This power law gives a reasonable fit although not as good as the Zener model for a slow indentation, with parameters: t 0 =0.056s, B=288.3, 0.104 (cf. blue curve in Fig. 8). In the case of a rapid indentation, the comparison between the Maxwell model with two relaxation times and the power law model is shown in Fig. 9. The obtained parameters were t 0 =0.016, =0.265, B=199; clearly the agreement with the experimental curve is less satisfactory for the power law model than for the Maxwell model with two relaxation times, especially at short times where the decrease of the force is too slow (insert of Fig. 9). Nevertheless it is worth noting that the value of the coefficient  is not very different from the one =0.20 found by a different experimental method [START_REF] Balland | Power laws in microrheology experiments on living cells: Comparative analysis and modeling[END_REF].

(1-,t 0 ) becomes (1-) so:     1 cos 2             GB ( 41
)
The red curve clearly shows the two relaxation frequencies of the generalized Maxwell model, and a high frequency plateau, whereas the power law model shows a continuously increasing modulus. The addition of more than two relaxation times in the generalized Maxwell model is difficult to perform because of too many parameters which can give similar residue for different set of parameters; it can improve a little bit the fit and gives a smoother frequency dependency (blue line in Fig. 11) but in any event it remains true that both at low and high frequencies the behavior of the two models remains very different. At low frequency we find clearly a zero frequency modulus, G 0 corresponding to a constant position of the indenter after the relaxation sequence; this also confirmed by the elastic recovery of the indenter position at the end of the cycle described in Fig. 2 whereas in the power law model there is a continuous decrease of the modulus with the frequency. Actually it is possible to add a zero frequency modulus in Eq.36 to take into account the solid like behavior, but at high frequency we find a plateau with the Maxwell model contrary to the power law model. This plateau is, of course, inherent to the Maxwell model and is imposed by the shortest relaxation time; it could be that experimentally it is not possible to test higher frequencies because the indentation velocity is limited. Actually, as can be seen in Fig. 7, we are well able to detect this plateau for indentation velocity higher than 10m/s. Furthermore, as can be seen in the insert of Fig. 10, the experimental relaxation is faster than the one predicted by the power law model, which means that this model does not well capture the physics of the relaxation process in the high frequency range. This frequency cut-off could be, for instance, related to a typical unbinding time of crosslinks, which blocks the relaxation of the actin network for smaller times [START_REF] Broedersz | Cross-Link-Governed Dynamics of Biopolymer Networks[END_REF].

Conclusions

In the section 2 of this paper it was shown, with the help of a comparison with FEM results on a viscoelastic solid, that a generalized Hertz theory was very well adapted to describe, not only the relaxation part, but the whole indentation curve, obtained with a spherical probe mounted on the tip of an AFM. Furthermore this equation takes into account the finite thickness of the cell with respect to the indentation depth. In the case of a Zener model, the differential equation relating the indentation to the vertical motion, z(t) of the piezoelectric transducer can be written explicitly: Eq.27 whatever the shape of z(t). For any creep function G(t), the integral equation Eq.23 can be used to obtain the parameters included in the function G(t). A proper determination of the contact between the probe and the cell and also the correction due to the viscous dissipation on the cantilever was proposed and described in Fig. 2. Based on these methods, the viscoelastic properties of Hep G2 cancer cells were obtained. It appears that high velocities of indentation reveal a short relaxation time which is generally ignored due to the use of smaller indentation velocity, This is for instance the case [2][34],

where the typical time of compression was several seconds which did not allow to observe this fast relaxation time. A generalized Maxwell solid with two relaxation times is shown to rather well describe the indentation function except for the very first part of the relaxation (t<0.1s). A power law function for G(t) gives a poorer agreement in this domain and furthermore is not able to represent the solid behavior at low frequency. A tentative interpretation of this result in terms of the actin network was proposed even if it was beyond the scope of this paper which was to give a firm basis for the determination of viscoelastic functions with the help of an AFM. These results will also be used in a future work to predict the motion of magnetic nanoparticles actuated at the surface of the cell by a magnetic field.

Fig. 1 :

 1 Fig. 1: The organization of actin (left) and microtubules (right) of Hep G2 shown by fluorescent microscopy

Fig. 2 :

 2 Fig. 2: Left is an indentation cycle; right is the method for the determination of the contact point and of the force/indentation curve. Upper graph: experimental force versus time. Lower graph: motion of the cantilever head versus time

Fig. 4 :

 4 Fig. 4: Comparison between FEM result and the solution of Eq. 27 for the indentation versus time in disk of different thicknesses.

Fig. 5 :

 5 Fig. 5 : Maximum and average stress with respect to the indentation depth for a purely elastic medium; the diameter of the bead is 5µm and the thickness of the layer is 10µm; the insert is the ratio of maximum stress to average stress

Fig. 6 :

 6 Fig. 6: AFM topography of Hep G2 cell, the point A to G are the measurements where has been carried out

Fig. 7 :

 7 Fig. 7: Young's modulus versus indentation velocity. In the insert: force versus indentation for two different velocities (v=100µm/s and 100nm/s) and their fit by Eq. 20

Fig. 9 :

 9 Fig. 9 : Rapid indentation: applied force versus time during an indentation with a

Fig. 10 :Fig. 11 :

 1011 Fig. 10 : Same conditions as in Fig. 7. Comparison between a fit with two Maxwell branches and a fit with a power law

Viscoelastic parameters deduced from a fit of G(t) with the generalized Maxwell model (Eq.29) with one Maxwell branch for slow indentation and two Maxwell branches for rapid indentation; the thickness of the cell was L=10m; the values between parenthesis are obtained for an infinite medium

  Different branches of the Maxwell model can tentatively be associated with different deformations of the actin network; the first spring with the shear modulus G 0 gives the solid like behavior and is likely mainly due to the enthalpic bending stiffness of the actin network.

		0	G 1	 1 
	Average	210 (247.6) 118.4 (128.3) 5.6 (6.1)
	Standard deviation 22.8 (28.5)	19.8 (20.7) 2.1 (2.2)
	Rapid indentation	G 0	G 1	 1 	G2	 2 
	Average	184.3 (214.5) 149.4 (168.6) 3.07 (3.84) 451.4 (538.5) 0.057 (0.065)
	Standard deviation 65.8 (74.0)	30.2 (35.1) 0.99 (1.16) 110.2 (142.7) 0.003 (0.018)
	Tab 1:

Acknowledgement

The authors thank the PACA Region and CNRS Cote d'Azur, for their financial support of the project BIOMAG and we are grateful to M. Girard from CEMEF (Sophia-Antipolis) for her help with Abaqus simulations and G.Bossis from IGMM (Montpellier) for providing the images with fluorescent microscopy.