Amira Barhoumi
email: amira.barhoumi.etu@univ-lemans.fr

Chafik Aloulou
email: chafik.aloulou@fsegs.rnu.tn

Nathalie Camelin
email: nathalie.camelin@univ-lemans.fr

Yannick Estève
email: yannick.esteve@univ-lemans.fr

Hadrich Belguith
email: l.belguith@fsegs.rnu.tn

Arabic Sentiment analysis : an empirical study of machine translation's impact

Keywords: Sentiment analysis, opinion mining, machine translation, document embeddings, machine learning, Arabic language

The largest amount of Sentiment Analysis has been carried out for English language. To deal with Arabic sentiment analysis, machine translation of English resources or Arabic texts may be applied to built Arabic sentiment analysis systems. In this paper, we translate Arabic dataset into English and study the impact of machine translation while considering a standard Arabic system as a baseline. Experiments show that sentiment analysis of Arabic content translated into English reach a competitive performance with respect to standard sentiment analysis of Arabic texts. This suggests that machine translation can successfully transfer the expression of sentiment or polarity. Moreover, we explored the multi-domain extending of training data in order to enhance performance and we show that we should have, in the training set, data whose domain is the same as the domain of evaluation dataset.

Introduction

Sentiment analysis (SA) is a classification task that involves building systems recognizing the opinion expressed in natural language sentences. Its goal is commonly to identify the subjectivity (objective/subjective) and the polarity (positive/negative) of a given text [START_REF] Pang | Opinion mining and sentiment analysis[END_REF]. Several research in SA have been carried out for English language. Nevertheless, there are few works that have been done for Arabic. This could be explained by the low number of resources developed in Arabic and their unavailability [START_REF] Al-Kabi | A prototype for a standard arabic sentiment analysis corpus[END_REF]. In sentiment analysis, making reliable resources is expensive and time-consuming because it needs experts to annotate corpora and built knowledge documents (lexicons, sentiment ontology, etc.). To avoid the cost of setting up resources, one track one may adopt is to translate existing resources from low-resources language to high-resources language and apply efficient methods for SA in that language.

Indeed, advances in machine translation (MT) have given Natural Language Processing (NLP) systems a strong boost. It was incorporated in many NLP applications such as online translation services, information extraction, document retrieval, etc. Research in opinion analysis took also advantage from machine translation, especially with under-resourced languages. In this paper, we focus on the task of Arabic Sentiment Analysis (ASA) studying the impact of machine translation from Arabic to English on sentiment texts. Throughout a set of experiments presented in this paper, we answer several research questions such as:

1. Do MT systems alter the subjectivity or polarity expressed in the source Arabic text? 2. What kind of performances can be reached using MT systems in SA compared to standard SA on Arabic texts? 3. What kind of interest is there to use MT in ASA? 4. Does the use of a larger training data improve the results even if additional data comes from an other domain?

In this paper, we investigate how sentiment is preserved after machine translation. We conducted several experiments to study: the impact of MT systems; and the use of additional data from other domain in the training set. The remaining of this paper is organized as follow. In the next section, we give a short overview of ASA systems with a special attention to those using MT systems. Section 3 describes the methodology we propose to measure the impact of the use of MT in ASA. Section 4 presents the experimental setup, evaluation, performance and results discussion. Finally, we conclude and make suggestions for future researches in section 5 .

Related works

Research on sentiment analysis show a great interest from the scientific community as showed by various evaluation campaign [START_REF]Proceedings of the 9th International Workshop on Semantic Evaluation (SemEval 2015). Association for Computational Linguistics[END_REF], [START_REF] May | Semeval-2016 task 8: Meaning representation parsing[END_REF], [START_REF] Rosenthal | SemEval-2017 task 4: Sentiment analysis in Twitter[END_REF], [START_REF] Mohammad | Semeval-2018 task 1: Affect in tweets[END_REF], [START_REF] Benamara | Analyse d'opinion et langage figuratif dans des tweets présentation et résultats du défi fouille de textes DEFT2017[END_REF]. If wide research has been carried out for English language, few works has been done for Arabic.

In this work, we focus on Arabic language. Indeed, many researchers have investigated sentiment analysis and opinion mining from different classification approaches. However, limited research is conducted on Arabic sentiment analysis. The research field on Arabic3 is characterized by a luck on sentiment resources: annotated corpora and lexicons 4 . Thus, many research on ASA use machine translation to built Arabic resources thanks to the recent progress of MT systems especially from other languages into English.

In this section, we present some works that use machine translation to build Arabic Sentiment Analysis systems. In this framework, two tracks are possible: either translate from Arabic to English and train a system on English texts, or translate English resources into Arabic and train a system on Arabic texts. [START_REF] Salameh | Sentiment after translation: A casestudy on arabic social media posts[END_REF] and [START_REF] Mohammad | How translation alters sentiment[END_REF] tested the two tracks and they improve performance in sentiment prediction with both manual and machine translation documents. They conclude that sentiment analysis systems are able to capture polarity from translated texts. Moreover, [START_REF] Elnagar | Comparative study of sentiment classification for automated translated latin reviews into arabic[END_REF] show that sentiment analysis of both original English dataset and the Arabic translated one produce comparable results. They infer that automatically translated English datasets shall be used as resources for building robust Arabic sentiment analysis systems. It's commonly known that building SA system requires corpora and/or lexicon resources. Some researches ( [START_REF] Refaee | Benchmarking machine translated sentiment analysis for arabic tweets[END_REF], [START_REF] Mohammad | How translation alters sentiment[END_REF]) translate Arabic corpora into English and others translate English lexicon into Arabic. In fact, [START_REF] Eskander | Slsa: A sentiment lexicon for standard arabic[END_REF] built an Arabic lexicon SLSA based on the translation of the English lexicon SentiWordNet5 and Arabic morphological analyzer. In the same way, [START_REF] Alotaibi | Extending the knowledge of the arabic sentiment classification using a foreign external lexical source[END_REF], [START_REF] Duwairi | Detecting sentiment embedded in arabic social media-a lexicon-based approach[END_REF] and [START_REF] Refaee | Benchmarking machine translated sentiment analysis for arabic tweets[END_REF] rely on the best state-of-the-art MT system and propose a method that leverages the available resource in another language such as the English SentiWordNet. Moreover, [START_REF] Guellil | Sentialg: Automated corpus annotation for algerian sentiment analysis[END_REF] exploited the machine translation to construct an Algerian sentiment lexicon (containing words in both Arabic and Arabizi). We can also mention [START_REF] Mahyoub | Building an arabic sentiment lexicon using semi-supervised learning[END_REF] that used a seed list contained 14 English words , translated them to Arabic and filtered them based on the Arabic WordNet. [START_REF] Ibrahim | Automatic expandable large-scale sentiment lexicon of modern standard arabic and colloquial[END_REF] do the same by adding synonyms and translation to a basic lexicon. Moreover, [START_REF] Al-Horaibi | Sentiment analysis of arabic tweets using semantic resources[END_REF] use both Arabic and English sentiment lexicons to classify the Arabic tweets into three sentiment categories (positive or negative or neutral).

The majority of works concludes that translation brings competitive results. However, [START_REF] El-Beltagy | Weightednileulex: A scored arabic sentiment lexicon for improved sentiment analysis[END_REF] has shown that sentiment analysis accuracy does suffer when using translated lexicon and the quality of such lexicon is not as high as a manually constructed one.

[4] explores cross-lingual sentiment classification from English to Arabic, without any manual annotation effort, and found that it is easy to build and does not require deep linguistic analysis. They conclude that a good classification model can be obtained from translated corpora regardless of the noise added by machine translation.

Methodology

In this work, we are interested in studying the impact of machine translation on Arabic sentiment analysis. Thus, we propose the following method: to automatically translate Arabic text into English, to learn a classifier based on that translated text using English document embeddings as input. As a consequence, we can explore how English translation of Arabic text alters or not the expression and/or detection of sentiments.

In an other hand, we consider a standard ASA system trained on Arabic texts as a baseline. We will compare both systems performances.

The whole experimental scheme is outlined below:

-Identify the arabic dataset dataset ar.

-Learn an SA system on dataset ar train based on pre-learned Arabic document embeddings. This system is denoted as the ASA-baseline. -Run the ASA-baseline system on dataset ar test and compute its performance. -Translate dataset ar from Arabic to English and obtain the English version dataset en.

-Learn an SA system on dataset en train based on pre-learned English document embeddings. This system is denoted ASA-MT. -Run the ASA-MT classifier on on dataset en test and compute its performance. -Compare the performance of both systems ASA-baseline on dataset ar and ASA-MT on dataset en and draw some inferences.

Experiments and results

In this work, we choose to focus the sentiment analysis only considering the polarity expressed. As a consequence, the two SA systems implemented are binary classifiers that predict the two classes: positive and negative.

Architecture

We explore two classifiers: logistic regression (LR) and multi-layer perceptron (MLP) 6 . The input vector of each classifier is the embedding obtained by learning paragraph vector algorithm. It allows obtaining distributed representations (Doc2vec) for any length sequence, ranging from phrases to documents. It efficiently computes document vector representations in a dimensional vector space. The Doc2vec representations were used for English sentiment analysis by Le and Mikolov [START_REF] Le | Distributed representations of sentences and documents[END_REF] who achieve the best performance with paragraph vector compared to other approaches on IMDB [START_REF] Maas | Learning word vectors for sentiment analysis[END_REF] dataset which contains 100000 film reviews.

Motivated by their work, we choose to use Doc2vec embeddings as input of our two systems. The input vector is a concatenation of two vectors: one learned from distributed memory version (DM) and one learned from distributed bag of words version (DBOW). Each one have 400 dimensions. So that, 800 is the dimension of the classifier input. As a result, we kept the same neural architecture and the same hyperparameters of paragraph vector model used in [START_REF] Le | Distributed representations of sentences and documents[END_REF].

In order to translate Arabic texts into English, we used the LIUM 7 machine translation system. It is based on the statistical machine translation engine Moses8 [START_REF] Koehn | Moses: Open source toolkit for statistical machine translation[END_REF]. It takes large quantities of parallel data (Arabic/English) and uses coocurrences of words and phrases to infer translation correspondences between the two languages. For decoding, Moses finds the highest scoring sentence in the target language (here, English).

Training data

The learning of Doc2vec representations needs a big corpus. According to our knowledge, LABR dataset [START_REF] Nabil | Labr: A large scale arabic sentiment analysis benchmark[END_REF] is the biggest arabic dataset for SA that is freely available 9 .

Hence, we used LABR corpus for a set of various ASA experiments. This corpus consists of 63257 book reviews written in modern standard Arabic (MSA) and colloquial Arabic. Each review is associated by its author with a rate ranging from 1 to 5 stars. Table 1 describes the distribution of the reviews among the different ratings. It is customary to consider reviews with 1 or 2 stars as negative reviews and those with 4 or 5 stars as positive. Reviews with 3 stars are neutral and they are not considered in classification. Thus, final corpus is reduced to 40845 reviews (68% positive) for the training corpus and 10211 for the test corpus (69% positive). Note that 10% of the training set is used as a development corpus.

Feature extraction and experimental setup

As mentioned in the architecture sub-section, two classifiers were investigated: logistic regression and multi-layer perceptron. The input of each classifier is a set of document embeddings.

We tested three different types of document embedding: DM, DBOW and the concatenation DM+DBOW. We trained the classifiers with different learning rates (10 -3 , 10 -4 et 10 -5 ). While varying different hyper-parameters, LR classifier gives better results than MLP. That is why we only reported in this work the performances of logistic regression, MLP results do not inferred more interesting information.

To evaluate the performance of SA on the LABR dataset, we carried out several experiments using various configurations. All the experiments were conducted in Python using Theano10 for classification and Gensim 11 for learning vector representation. The performance is measured with error rate 12 metric which calculates the percentage of misclassified examples.

Results and discussion

First, the two system ASA-baseline and ASA-MT were running as detailed in previous sub-sections. Table 2 shows the error rates of both systems obtained with the logistic regression classifier. The error rate of ASA-baseline system is 25.37% that drops to 23.7% when using the ASA-MT system. So, we note a gain of 1.7% with the proposed ASA-MT system. This latter is trained on LABR en dataset which represent the English translation of LABR ar used in the ASA-baseline. So that, we highlight that machine translation does not seem to alter the expression of a polarity and brings a competitive (even better) accuracy with respect to Arabic sentiment analysis.

In order to explain and understand why translating the arabic LABR dataset in English enhance the performance, we analyzed the translated corpus LABR en. One transformation done during the translation process is to remove all nontranslated words from LABR en as these words may be ambiguous or noisy in the English dataset. Considering that choice and the observed results, we have made the assumption that non-translated words may carry ambiguity. To validate this hypothesis, we conducted a non-translated Arabic word's analysis. It shows that these non-translated words are mainly: dialectal words, proper names, typo words, words with repetition of characters, non-Arabic origin words written with Arabic letters (such as: /Alrfyw/ for review, /Albrwtk$n/for protection, etc.). Following this analysis, we consider that these non-translated words seem to disrupt the polarity detection, so we choose to consider these words as noisy ones. In order to see if sufficient information is kept in the English corpus to extract the polarity, we choose to remove the noisy words from the Arabic dataset. This corpus is denoted LABR ar MTremove. A new LR classifier is learned on LABR ar MTremove. We hope that this new classifier would have a positive impact on performance. Unfortunately, we obtained an error rate of 26.86% in that configuration which is greater than 25.37% initially obtained on LABR ar. In that way, this removal words seems to be informative for the standard ASA-baseline system, even if they are not for the MT-system.

Other ASA systems using NLP preprocess tools were implemented like [START_REF] Barhoumi | Document embeddings for arabic sentiment analysis[END_REF] that chose to apply light stemming as a pre-processing step on LABR ar. As a result, they obtained an error rate of 23.31% which is better than 25.37% obtained by our ASA-baseline. This means that light stemming is a reliable pre-process step for Arabic sentiment analysis. Moreover, it is higher (approximatly equal) than 23.7% obtained on the English dataset LABR en. We infer that machine translation, as a statistic tool, and light stemming, as a linguistic tool, almost behave the same way for the Arabic sentiment analysis task. So for languages without such linguistic tools, we could apply machine translation to have a good sentiment analysis system, especially with the progress of machine translation field. In other words, we could use machine translation if linguistic specific tool is not yet developed or not powerful.

To go further to explore the impact of MT systems in ASA, we believe that the performance obtained with the proposed ASA-MT system may be enhanced using a larger amount of training data. Moreover, we believe that a larger corpus of reviews, whatever the domain, should be positive. In this perspective, we made some experiments while varying the domain (books, films) of reviews. This setup is as follows:

-Training set composed of IMDB dataset [START_REF] Maas | Learning word vectors for sentiment analysis[END_REF]

Conclusion and future works

In this paper, we presented a set of experiments to study the impact of English machine translation on sentiment analysis of Arabic reviews. Our experiments show that sentiment analysis of English translation is better than sentiment analysis of gross Arabic texts. We observed that machine translation does not alter the polarity prediction. Moreover, we found that sentiment analysis of English translation reach competitive results with respect to sentiment analysis of light-stemmed Arabic texts. So, we could generalize, whatever the language, that machine translation could be used if such linguistic tools (light stemming) do not exist or are not efficient. We have also explored the track of extending the training corpus and we proved the interest of keeping data whose domain is the same as the test set domain. We also showed that mixing the training data domains performs better than simply changing the training domain.

As future work, we think combining light stemming and machine translation based systems. In fact, these systems behave differently as they use different techniques. We strongly think that combining them could enhance performances.

Moreover, one other perspective is to investigate the use of deep learning classifiers, especially convolutional neural networks (CNN) which are efficient in English sentiment analysis ( [START_REF] Kim | Convolutional neural networks for sentence classification[END_REF], [START_REF] Alayba | A combined cnn and lstm model for arabic sentiment analysis[END_REF], [START_REF] Zhang | Dependency sensitive convolutional neural networks for modeling sentences and documents[END_REF], [START_REF] Croce | Injecting sentiment information in contextaware convolutional neural networks[END_REF]). We think testing CNN and other deep learning techniques on ASA.

Finally, we would explore different types of embedding, in addition to document embeddings, like sentiment-specific word embeddings [START_REF] Yu | Refining word embeddings for sentiment analysis[END_REF].

Table 1 .

 1 Distribution of LABR dataset among rating stars.

		Very negative Negative Neutral	Positive	Very positive Total
	Training 2331	4195	9762	15189	19129	50606
	Test	608	1090	2439	3865	1649	12651

Table 2 .

 2 Error rates of Logistic regression classifier over Arabic and English-translated datasets: LABR ar and LABR en respectively.

	System	Dataset input	Error rate
	ASA-baseline	LABR ar	25.37%
	ASA-MT	LABR en	23.70%

  . -Training set composed of IMDB dataset and the LABR training set. -Training set composed of LABR training set and the classifier's parameters are initialized with the parameters obtained after training the classifier on IMDB dataset.Results are reported in Table3. We can observe that, by changing the domain of the training set, the error rate increases from 23.7% to 29.33%. However, while adding reviews of other domain to the training set, we obtain 26.94% as error rate. So, it performs better than changing completely the training domain (26.94% vs. 29.33%). Other than mixing data with different domains, multi-domain experiment can be also established by initializing the classifier parameters with those obtained while training on dataset of other domain. The error rate reaches 24.84% with this technique. It is better than training on domain-mixed data (24.84% vs. 26.94%).

Table 3 .

 3 Experiment's configurations to enlarge the training set.

		Partition sets	Error rate
		Train = IMDB
	Domain change	Dev = LABR-dev	29.33%
		Test = LABR-test
		Train = IMDB + LABR-train
		Dev = LABR-dev	26.94%
	Multi-Domain	Test = LABR-test Train = LABR-train
		Dev = LABR-dev	24.84%
		Test = LABR-test

For an overview of Arabic sentiment analysis field,[START_REF] Al-Ayyoub | A comprehensive survey of arabic sentiment analysis[END_REF] build a survey.

[START_REF] Al-Kabi | A prototype for a standard arabic sentiment analysis corpus[END_REF] and[START_REF] Boudad | Sentiment analysis in arabic: A review of the literature[END_REF] summarizes all freely available corpora for Arabic sentiment analysis task

http://sentiwordnet.isti.cnr.it/

The MLP contains 3 layers: the input layer whose number of neurons is equal to the size of the input vector, one hidden layer with 50 units and the output layer with 2 neurons to predict the polarity of the input text: either positive or negative

https://lium.univ-lemans.fr/

http://www.statmt.org/moses/

LARB dataset is available on http://www.mohamedaly.info/datasets/labr

http://deeplearning.net/software/theano/

https://radimrehurek.com/gensim/


Cariou S 1°

Ricoux P 4°

Assumani B 6°

Reliability Improvement of Odour Detection Thresholds Bibliographic Data

Keywords: Odour Detection Thresholds (ODT), Data mining, Reliability, Completeness, Uncertainty

Odour control is an important industrial and economical issue as it is a criterion in purchase and use of a material. The minimal concentration of a pure compound allowing to perceive its odour, called Odour Detection Threshold (ODT), is a key parameter of the odour control. Each compound has its own ODT. Literature is the main source to obtain ODT. Nevertheless, there are a lot of compounds with no reported ODTs and when ODTs are available, they are marred by a high variability. Another expensive and time-consuming way to obtain ODT is the measurement. This paper proposes a validated cleaning methodology to reduce uncertainty of available ODTs. This methodology will be consolidated by our own experimental measurements. Next, we predict missing ODTs as a function of chemical and physical variables. The proposed cleaning methodology leads to eliminate 39% of compounds with at least one ODT while conducting 84% of positive scenarios (on 37 compounds). The missing ODTs are predicted with an error of 0.83 for the train and 1.14 for the test (on a log10 scale). Given the uncertainty of data, the model is sufficient. This approach allows working with a lower uncertainty on available ODTs and predicts missing ODTs with a satisfactory model.

Introduction

In the industrial environment, there is a growing need to identify compounds responsible for an unpleasant odour. This identification depends firstly on the Odour Detection Threshold (ODT) of each compound. We define ODT as the minimal concentration of a pure compound allowing to perceive its odour. The principal source of ODTs is the literature [START_REF] Leonardos | Odor Threshold Determinations of 53 Odorant Chemicals[END_REF][START_REF] Ruth | Odor Thresholds and Irritation Levels of Several Chemical Substances : A Review[END_REF][START_REF] Devos | Standardized human olfactory thresholds[END_REF][START_REF]ASTM: Odor Thresholds for Chemicals with Established Occupational Health Standards[END_REF][START_REF]Reference guide to odor thresholds for hazardous air pollutants listed in the clean air act amendements of 1990[END_REF][START_REF] Van Gemert | Odour Thresholds[END_REF]. But, in literature, there are a lot of compounds with no reported ODTs and when ODTs data are available, they present a high variability. This situation implies a high uncertainty of ODTs. This variability can be illustrated with the butyl acetate example whose ODT values range from 0.030 mg/m 3 to 480 mg/m 3 for 14 publications [START_REF] Van Gemert | Odour Thresholds[END_REF]. The variability could potentially be explained by a set of parameters such as: difference of methods; existence and year of normalization; sample quality; environmental conditions, culture of authors; panel selection and their intrinsic diversity; panel correction, etc. Given all these potential sources of variability, there is a need to improve the reliability of these data. Unfortunately, the sparsity of these data makes most of the potential interesting statistical tools unusable. This sparsity comes from non-uniformity of the information from one author to another. Several papers have already highlighted this issue [START_REF] Ruth | Odor Thresholds and Irritation Levels of Several Chemical Substances : A Review[END_REF][START_REF] Devos | Standardized human olfactory thresholds[END_REF] and have tried to compare publications [START_REF] Patte | Selected and standardized values of suprathreshold odor intensities for 110 substances[END_REF] but usually, researchers make a subjective sorting or simply use the mean of values [START_REF] Devos | Standardized human olfactory thresholds[END_REF][START_REF]Reference guide to odor thresholds for hazardous air pollutants listed in the clean air act amendements of 1990[END_REF]. But, in any event, even if dataset is used without cleaning, neither the mean or the geometric mean nor the density application is justified if the ODT values distribution is not identified. As a consequence, our paper aims to find relevant methodology allowing improving the reliability of these data. This improvement gathers an approach aiming to decrease the variability of available ODTs by a cleaning methodology. Next, we complete not reported compounds by a predictive modelling of ODTs as a function of chemical / physical variables. Explicative variables generated by the methodology are precise opposed to bibliographic ODTs.

Methodology

Software

The ODT values from the literature are collected in an Excel sheet and the reliability improvement is realized on the R software.

Database Construction

To analyze ODT, a database is constructed (Fig. 1). This database gathers quite a few volatile and odorous compounds. These compounds constitute the rows of the database.

For each compounds, a state of art of ODTs is done to identify publications containing ODT information [START_REF] Van Gemert | Odour Thresholds[END_REF]. These publications constitute the columns. This database will be the support of the uncertainty reduction and prediction. To date, it contains 161 publications (columns) and 1 254 compounds (rows) including only 650 compounds with at least one ODT available.

First Step of Cleaning: Provisional Publications' Isolation

The first step consists in eliminating the least reliable studies. Ideally, it would have been relevant to consider only studies containing a reliable repeatability of measurements. In this way the Cochran test based on the standard deviation of each article would have been applied [11]. Unfortunately, this information is rarely available. Considering only studies with available repeatability is therefore too restrictive. Consequently, we use the criterion of number of ODT measurements realized per study. Indeed, we have noticed that the higher the number of ODT determinations in a publication, the better the technique of analysis was described, and hence the more reliable the results. It was decided to eliminate ten percent of the total values of the database and that corresponds to publications containing less than four ODT values (Fig. 2). To define an outlier in this context, the ODT values distribution has to be defined. On the basis of an expert observation on well-known molecules, the hypothesis is that for each given compound, ODTs follow a normal distribution.

Considering the low amount of ODT values per compound (Fig. 3), the most suitable solutions to observe this distribution, are the Kernel Density Estimation (KDE) and the QQplot representation of the eleven compounds with more than ten ODT values after the first step of the cleaning.

Fig. 3. Reported ODT frequency

The KDE is calculated with the "geom_density" function and the QQplot with "qqnorm" function. The correlation coefficient (CC) of QQplot of these compounds is calculated (Fig. 4). The inter-laboratory reproducibility Grubbs test is applied [11]. The confidence level chosen is 95%.

Outliers Values Elimination According to Normal Distribution

The Grubbs test is not applicable with less than three values [START_REF]ASTM International: ASTM E178-00 Standard practice for dealing with outlying observation[END_REF]. In this way, only compounds measured more than twice can be compared with the rest of the database. Thereafter, ODTs of these compounds are considered as the "tested" ODTs. That implies that some compounds' relevance cannot be analyzed at this step.

Another rule of comparison has to be set up. This is an extension of the cleaning methodology to the elimination of publications outliers. This next rule was guided by the fact that, at the previous step, the eliminated ODT values often belong to the same publications.

Third Step of Cleaning Methodology: Publication Outliers Elimination

At this step, the relevance of each publication is measured by the number of exclusion of that article at the second step. An index is calculated for each author: the exclusion frequency ratio in equation [START_REF] Leonardos | Odor Threshold Determinations of 53 Odorant Chemicals[END_REF]. The ODT values of publications have not been compared in the same way. That's why a "tested rate" is calculated. For each publication, this "tested rate" is defined as the percentage of tested ODTs among all of the ODTs of this publication. Indeed, the higher the "tested rate" is; the better the 𝐸𝐹's reliability is.

First, publications with "tested rate" under 1/3 are eliminated. Then, after a Hierarchical Ascendant Classification (HAC) applied on the EF of selected publications, the Ward's distance index [START_REF] Ward | Hierarchical grouping to optimize an objective function[END_REF] is used. This classification allows to statistically separate publications based on EF values. The dendrogram and the SPRSQ graph are used to determine the number of groups. The group containing the lowest EF is finally retained.

Validation of the Cleaning

ODT of 44 compounds were measured experimentally in our laboratory to validate the relevance of this cleaning. These ODT were measured using the norm EN 13725 [14]. For these 44 compounds, whose ODT values have been measured experimentally, 40 were reported in the literature. For these 40 compounds, three barycenters are defined:

the one of raw ODT values of the literature (X1), the one of remaining ODT values after the cleaning methodology application (X2) and the one of our experimental ODT values (X3). And then two differences are calculated: X1-X3 (Z1) and X2-X3 (Z2). These two differences are compared. This comparison allows to observe if the cleaning methodology leads to approach the experimental result.

To state on the relevance of the cleaning the criterion to select the "ideal" case is the mean of the log10 of the confidence interval obtained thanks the repeatability level of our experimental measurements.

Completeness of the Database by Predictive Modelling

The completeness of the database is based on a Quantitative Structure-Property Relationship (QSPR) approach. This approach consists in predicting a variable (Y) as a function of chemical and physical certain variables (X). Here the Y variable is the mean value of log10 values of cleaned ODTs. Explicative variables X are essentially calculated from the structure of compounds (66 variables). There are compositional and constitutional indexes (the number and nature of atoms, the molecular weight, the unsaturation), topological indexes which the majority are defined in Todeschini's publication [START_REF] Todeschini | Molecular Descriptors for Chemoinformatics[END_REF] and electrotopological indexes as ZEP index [START_REF] Berinde | On a matrix representation of molecular[END_REF]. The model approach is then divided in 2 principal steps. The first is a reduction of the number of significative input variables by a lasso technique. To predict Y as a function of X, a classical Support Vector Machine (SVM) model was constructed. Some "bad" observations was eliminated from this learning process.The partition of the dataset between the train and test sets is made in order that the train represents 75% of the dataset and reflects the variance of this one.

Results and Discussion

Visualization of the Proportion of Missing Values of the Database

Missing value proportion of the database (Fig. 1) is presented as a heatmap (Fig. 5). The total proportion of missing values is 99%. This proportion enforces the inability to use statistical classification tools and implies the necessity to implement a more reliable methodology.

Result of the First Step

For the three "data cleaning" steps, information loss is presented from three perspectives (Table 1): the number of publications, the number of ODT and the number of compounds with at least one ODT. This choice is made because the most important information to monitor is the number of compounds with at least one ODT. Indeed, as it was presented, we want to predict the ODT behavior of all the compounds of the database (1254) with reported ones in literature. That's why, it's important to keep a satisfactory proportion of reported compounds. The first step leads to isolate 55% of the publications. Even if half of the publications have been eliminated, only 3% of compounds with at least one ODT have been eliminated. This proportion is acceptable. The lack of reliability on these sources justifies to reject them.

Result of the Second Step

The Grubbs test is applied on compounds which ODT have been measured at least 3 times. Results of the information loss of this step are presented in Table 2. There is only a slight loss of information on the ODT values (10%) and the two other rates are still constant. Furthermore, among this proportion, there is an average of 1.65 values eliminated per compound (over 80 compounds). This low percentage supports the use of a normal distribution.

Result of the Third Step

A summary of ODT elimination at the second step on all the 72 publications, will allow applying the third step. The focus is made on the Exclusion Frequency (EF) defined by the equation ( 1) and the "tested rate" defined in the 2.5 section.

As it was mentioned in this section, publications with "tested rate" under 1/3 are eliminated. A first information loss summary is done after the elimination of these publications with a "tested rate" lower than 33% (Table 3). The dendrogram of the HAC classification is applied on the EF of publications with a "tested rate" higher than 33% (Fig. 6).

It can be noticed that a high percentage of publications is eliminated (79%) but it represents almost the double of compounds with at least one ODT value in the literature (39%). The amount of total eliminated data is 53%. These rejected data have been considered irrelevant by our specific approach.

Validation of the Cleaning Methodology with Measured ODT

After the cleaning, the ODT value of 3 compounds out of 40 was totally eliminated and we make the comparison on the 37 remaining compounds. The mean of the repeatability of our measurement calculated in a log10 scale is 0.4. The ODTs obtained are then compared to this value as it was explained in the 2.6 section (Fig. 7). Many scenarios are recorded (Table 6) thanks to the graph of comparisons (Fig. 7). Initially, 16 compounds were in the range of the measured ODT. After the cleaning, there are 21 compounds (out of 37 compared compounds). This is a first positive added value of the cleaning methodology. Furthermore, it can be considered as positive situation after the cleaning the scenarios 1, 2, 3, 4 and 6 (Table 6). These cases represent 84% of positive cases (31 compounds out of 37). The worst scenario is the number 8. Despite of 6 compounds, this cleaning seems appropriate to decrease the uncertainty of available bibliographic ODT data.

ODTs Prediction

Applying lasso technique leads to eliminate 22 variables (out of 66) and SVM leads to a model with an error, in a log10 scale, of 0.83 on the train and 1.14 on the test with compounds eliminated (Fig. 8). As a matter of fact, among all the predictions, only less than 2.5% of observations was abnormally predicted in comparison with the others.

Fig. 8. ODTs predicted in function of ODTs observed

Given the uncertainty of ODT of most of compounds, the model is satisfactory as it predicts the ODT with a log10 error about 1. To our knowledge this is the first global model to predict ODT as a function of molecular characteristics on this quantity of compounds. This is an encouraging result which has to however be improved thereafter.

Conclusion

In conclusion, the proposed cleaning methodology eliminates 39% of compounds with at least one ODT while conducting to 84% of positive scenarios on ODT values (validation on 37 compounds). The missing ODTs are predicted with an error of 0.83 for the train and 1.14 for the test (on a log10 scale). Considering the data, it's a satisfactory model. This approach allows working with a low uncertainty on available ODTs and predicts missing ODTs with a satisfactory model. This is an encouraging methodology.

ODTs predicted in fucntion of ODTs observed

Thereafter, it would be interesting to strengthen the cleaning methodology with more measurements in our laboratory or a probabilistic validation. The global predicting model of ODT seems satisfactory to odour experts. It will be improved further by complementary local approach and / or expert knowledge

Fig. 1 .

 1 Fig. 1. ODT values (mg/m 3 ) database format

Fig. 2 . 1 .

 21 Fig. 2. ODT values loss rate in function of publications eliminated

Fig. 4 .

 4 Fig. 4. Correlation coefficient of the QQplot of compounds with more than 10 ODTOn this Fig.4the higher the correlation coefficient is, the more likely the distribution is normal. These results encourage applying a normal test on values to detect outliers.

  of the jth publication 𝑵 𝑶𝑫𝑻 𝒆𝒙𝒄𝒍𝒖𝒅𝒆𝒅 𝒋 : Number of ODT values of the jth publication excluded by Grubbs test 𝑵 𝑶𝑫𝑻 𝒕𝒆𝒔𝒕𝒆𝒅 𝒋 : Number of ODT values tested in the jth publication[START_REF] Leonardos | Odor Threshold Determinations of 53 Odorant Chemicals[END_REF]

Fig. 5 .

 5 Fig. 5. Proportion of missing ODT values in the literature (a red zone corresponds to information presence; a white zone corresponds to a lack of information)

Fig. 7 .

 7 Fig. 7. Comparison of the cleaned values and the raw values of literature with our experimental values (red: cleaned ODT difference; blue: raw ODT difference)

Table 1 .

 1 Information loss after the first step

			Cleaning at the first step
		Before	After	Information loss
	Nos. of publications	161	72	rate 55%
	Nos. of ODT	1501	1367	9%
	Nos. of compounds with at least 1 ODT	650	631	3%

Table 2 .

 2 Information loss after the second step by Grubbs test

			Cleaning at the second step
		Before	After	Information loss
	Nos. of publications	72	72	rate 0%
	Nos. of ODT values	1367	1235	10%
	Nos. of compounds with at least 1 ODT	631	631	0%

Table 3 .

 3 Information loss after the elimination of publicationss with tested rate lower than 33%

	Cleaning at the third step -Part 1
	Before	After	Information loss rate

Table 6 .

 6 Summary of the evolution by the cleaning methodology

	Scenario 1		… even if it was already correct	5 cases
	Scenario 2	Improvement of the value …	… which leads to a correct value 6 cases
	Scenario 3		… but the value still not correct	10 cases
	Scenario 4	No improvement of the value	… but it was already correct	7 cases
	Scenario 5	…	… and still not correct	3 cases
	Scenario 6		… but it still correct	3 cases
	Scenario 7	Degradation of the value …	… but it still not correct	1 cases
	Scenario 8		… becomes not correct	2 cases

Thanks to the dendrogram and the SPRSQ graph, 3 groups were realized (Fig. 6). The information loss summary is done after the second part of the third step (Table 4). This third step is rather drastic because it eliminates 37.7 % of compounds with at least one reported ODT value in the literature. Nevertheless, 393 reported compounds still remain, spread over 33 publications. The size of the database decreases with the increase of its reliability. Of course, we consider that it is better to work with less quantitative but more informative data.

Summary of the Cleaning Methodology

After the three cleaning methodology steps, the information loss statement is presented in Table 5.


Azadeh Sheidaei

Morteza Kazempour

Ali Hasanabadi

Fayyaz Nosouhi

Martine Pithioux

Daniel George
email: george@unistra.fr

Influence of bone microstructure distribution on developed mechanical energy for bone remodeling using a statistical reconstruction method

Keywords: Bone microstructure, statistical reconstruction, finite element model, mechanical behavior

The development of a predictive model for bone remodeling is becoming increasingly important for medical applications such as bone surgery or bone substitutes like prostheses. However, as bone remodeling is a complex multiphysics phenomenon and difficult to quantify experimentally, predictive numerical models remain, at best, phenomenologically driven. Patient dependency is often ignored as its influence is usually considered secondary, although it is known to play an important role over long periods of time. Another difficulty to study this patient dependency is the availability of experimental samples to carry out extensive analyses. Using our recently developed statistical reconstruction framework, a set of ''bone like'' microstructures with variety of distributions has been created to study pseudo ''patient variabilities. '' The method provides similar effective stiffness tensor, equivalent stresses, and strain energy distributions for the original and the statistically reconstructed samples. The main outcome of this study is the correlation of similar effective mechanical properties between samples when bone remodeling will depend on the local strain energy distribution as a function of each bone microstructure. It is expected that two different microstructures with equivalent bone volume fraction will lead to identical bone remodeling in a short period of time, whereas this needs to be proven for long term evolution. This work could be used to develop precise predictive numerical models while developing parametric studies on an infinite number of virtual samples and correlating patient dependency with more precise mechanobiological numerical models.

Introduction

It is well established in the literature that the interaction between mechanics and biology plays a crucial role in the interpretation of growing tissue behavior [START_REF] Carter | Mechanical factors in bone growth and development[END_REF][START_REF] Casanova | Temporal evolution of skeletal regenerated tissue: what can mechanical investigation add to biological?[END_REF][START_REF] Huiskes | Effects of mechanical forces on maintenance and adaptation of form in trabecular bone[END_REF]. This mechanobiology depends on the natural tendency of the biological tissue to develop optimizing behavior to gain proper mechanical strength with the minimum amount of material. For bone adaptation and remodeling, this optimization problem is at the origin of many studies since the early work by Wolff [START_REF] Wolff | Das Gesetz der Transformation der Knochen[END_REF] followed by many others [START_REF] Cowin | Wolff's law of trabecular architecture at remodeling equilibrium[END_REF][START_REF] Frost | Bone ''mass'' and the ''mechanostat'': a proposal[END_REF][START_REF] Lekszycki | Modeling of bone adaptation based on an optimal response hypothesis[END_REF][START_REF] Lekszycki | Functional adaptation of bone as an optimal control problem[END_REF][START_REF] Pivonka | Model structure and control of bone remodeling: a theoretical study[END_REF] exploring the internal architecture and external conformation of bones changes in accordance with mathematical laws. The evolution of bone microstructure being dependent on cells and microstructure distribution, various mechanical stimuli have been proposed to trigger bone adaptation [START_REF] Hegedus | Bone remodeling II: small strain adaptive elasticity[END_REF][START_REF] Andreau | Optimal-tuning PID control of adaptive materials for structural efficiency[END_REF][START_REF] Prendergast | Prediction of bone adaptation using damage accumulation[END_REF][START_REF] Doblare´ | Anisotropic bone remodelling model based on a continuum damage-repair theory[END_REF]. While continuum models, for the description of natural bone regeneration, are widely spread in the scientific literature, scale related effects of the material, are not often integrated in the structural evolution due to its highly heterogeneous nature at the small scales. At these scales, the classical continuum theory does not allow to describe the correct material behavior where the prediction made by Cauchy continuum theory starts being far away from the experimental evidence [START_REF] Dell'isola | At the origins and in the vanguard of peridynamics, non-local and highergradient continuum mechanics: an underestimated and still topical contribution of Gabrio Piola[END_REF][START_REF] Misra | Identification of higher-order elastic constants for grain assemblies based upon granular micromechanics[END_REF][START_REF] Placidi | Gedanken experiments for the determination of two-dimensional linear second gradient elasticity coefficients[END_REF][START_REF] Abali | Theory and computation of higher gradient elasticity theories based on action principles[END_REF][START_REF] Dell'isola | Higher gradient continua: the legacy of Piola, Mindlin, Sedov and Toupin and some future research perspectives[END_REF]. Hence, adequate homogenization methods [START_REF] Re´mond | Applied RVE reconstruction and homogenization of heterogeneous materials[END_REF][START_REF] Martin | A thermodynamically consistent model of bone rotary remodeling: a 2D study[END_REF] need to be developed in order to integrate all these local phenomena at the global scale.

Historically, bone mechanical behavior was dealt with using classical optimization processes based mainly on the developed mechanical strain energy density (SED) through the material [START_REF] Madeo | A continuum model for the bio-mechanical interactions between living tissue and bio-resorbable graft after bone reconstructive surgery[END_REF][START_REF] Lekszycki | A mixture model with evolving mass densities for describing synthesis and resorption phenomena in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Madeo | A second gradient continuum model accounting for some effects of microstructure on reconstructed bone remodeling[END_REF][START_REF] Andreaus | A 2-D continuum model of a mixture of bone tissue and bio-resorbable material for simulating mass density redistribution under load slowly variable in time[END_REF][START_REF] Scala | Mechanically-driven bone remodeling simulation: application to LIPUS treated rat calvarial defects[END_REF][START_REF] Giorgio | A visco-poroelastic model of functional adaptation in bones reconstructed with bio-resorbable materials[END_REF][START_REF] Spingarn | Multiphysics of bone remodeling: a 2D mesoscale activation simulation[END_REF][START_REF] Cuomo | Forms of the dissipation function for a class of viscoplastic models[END_REF][START_REF] Giorgio | Viscous second gradient porous materials for bones reconstructed with bioresorbable grafts[END_REF]. However, to be able to obtain a good predictive numerical model, the biological phenomena [START_REF] Burr | Basic and applied bone biology[END_REF] need to be linked with the surrounding physics [START_REF] Pivonka | Mathematical modeling in bone biology: From intracellular signaling to tissue mechanics[END_REF][START_REF] Lemaire | Multiphysical modelling of fluid transport through osto-articular media[END_REF][START_REF] Lemaire | What is the importance of multiphysical phenomena in bone remodelling signals expression? A multiscale perspective[END_REF][START_REF] Lemaire | Three-scale multiphysics modeling of transport phenomena within cortical bone[END_REF][START_REF] Sansalone | On the uncertainty propagation in multiscale modeling of cortical bone elasticity[END_REF], and account for the growing capillaries and nutrient supply [START_REF] Bednarczyk | E A novel mathematical model for growth of capillaries and nutrient supply with application to prediction of osteophyte onset[END_REF][START_REF] Lu | A novel coupled system of non-local integro-differential equations modelling Young's modulus evolution, nutrients' supply and consumption during bone fracture healing[END_REF][START_REF] Moya | Human mesenchymal stem cell failure to adapt to glucose shortage and rapidly use intracellular energy reserves through glycolysis explains poor cell survival after implantation[END_REF][START_REF] Paquet | Oxygen tension regulates human mesenchymal stem cell paracrine functions[END_REF]. Lately, a new form of mechanobiological stimulus for bone remodeling has emerged integrating coupled biological and mechanical effects accounting for cell migration and the nutrient supply chain [START_REF] George | Mechanobiological stimuli for bone remodeling: mechanical energy, cell nutriments and mobility[END_REF][START_REF] George | A multiphysics stimulus for continuum mechanics bone remodeling[END_REF][START_REF] Allena | Reaction-diffusion finite element model of lateral line primordium migration to explore cell leadership[END_REF]. Although this allows accounting for the couplings between the mechanics and biological phenomena, it remains a macroscopic approach and requires more insight about the locally occurring mechanobiology.

Bone microstructure distribution is optimized in order to sustain the externally applied mechanical loads and its evolution changes over the long term as a function of local biological information. However, on a macroscopic scale, where the bone density is characterized only as a homogeneous media, the influence of this microstructure over the long term remodeling kinematics is little investigated. Although some numerical models tried to predict the bone remodeling kinematics from known trabecular bone distributions [START_REF] Hollister | A homogenization sampling procedure for calculating trabecular bone effective stiffness and tissue level stress[END_REF][START_REF] Tsubota | Computer simulation of trabecular remodeling in human proximal femur using large-scale voxel FE models: approach to understanding Wolff's law[END_REF][START_REF] Jang | Computational simulation of simultaneous cortical and trabecular bone change in human proximal femur during bone remodeling[END_REF][START_REF] Marzaban | Numerical simulation of load-induced bone structural remodeling using stress-limit criterion[END_REF], these do not provide insight about long term evolutions. In addition, patient dependent bone microstructure is very variable [START_REF] Kersch | Measurement of structural anisotropy in femoral trabecular bone using clinicalresolution CT images[END_REF][START_REF] Lian | Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches[END_REF], this could impact directly the numerical model predictions. Lately, some numerical models have investigated the correlations between the bone stiffness and simplified microstructure distribution [START_REF] Goda | A micropolar anisotropic constitutive model of cancellous bone from discrete homogenization[END_REF][START_REF] Rath | Scaling relations between trabecular bone volume fraction and microstructure at different skeletal sites[END_REF][START_REF] Goda | A 3D elastic micropolar model of vertebral trabecular bone from lattice homogenization of the bone microstructure[END_REF][START_REF] Wierszycki | A two-scale approach for trabecular bone microstructure modeling based on computational homogenization procedure[END_REF][START_REF] Goda | Optimal internal architectures of femoral bone based on relaxation by homogenization and isotropic material design[END_REF] through multiscale homogenization procedures. However, the bone remodeling processes accounting for these microstructure variabilities at the macroscopic scale are not studied since these works are mostly related to the structural mechanics homogenization process.

The ultimate goal of this research is to study the influence of the local bone microstructure distribution on the macroscopic bone density evolution over long periods of time. The first step of this study require an adequate method for bone microstructure reconstruction in order to study large sets of specimens as it is not possible to extract that amount of real bone samples, and show that the proposed reconstruction method is able to provide corresponding results with different structural samples as equivalent effective stresses, strain energy and stiffness tensor. For this, we present here a methodology to study the influence of the bone microstructure distribution on the developed mechanical energy field. Real bone microstructures are extracted from computed tomography (CT) scan data and segmented/reconstructed to obtain the corresponding finite element (FE) numerical model. A statistical method is used to reconstruct several equivalent bone microstructures with identical homogeneous stiffness. FE models are then constructed from these statistical microstructures. The same mechanical loads are applied on the different FE models to extract the variabilities of the local SED that will be at the origin of the bone remodeling kinetics.

Experimental details

Specimen preparation

The real bone specimen was extracted from the femur of an 85-year-old woman post mortem. The adult bone was harvested from the corpse at the Timone Hospital (Marseille, France). The specimen was cut into a parallelepiped rectangle of length, width, and thickness dimensions 9.88 mm, 6.13 mm, and 2.29 mm, respectively, using a water-cooled low-speed diamond saw (Buehler Isomet 4000, Buehler, Lake Bluff, IL, USA). The length of each specimen was measured with a digital calliper (Absolute digimatik solar, Mitutoyo, Kanagawa, Japan, measurement error of 0.03 mm). The specimen was stored under vacuum at 220°C.

Image acquisition

The specimen was imaged at 5.5 mm voxel size using a micro-CT system (Phoenix Nanotom 180 S, General Electric, Germany). The scanning was done with a hydrated bone specimen, immersed in salt water inside a 25 mm inner diameter plastic tube. The axial specimen main axis (i.e., Haversian canal's principal orientation and medullar canal orientation) was aligned with the rotational sample holder axis. Scans were done with a field of view of 2284 × 2304 pixels, with a maximum voltage source of 90 kV, a current of 120 mA and a rotation of 360°providing 2400 projections per scan. The exposure time was 500 ms with a total scan time of 85 minutes. The bone density was measured at 92% of the overall specimen volume.

Theoretical modeling

The statistical reconstruction approach enables to obtain various different microstructures with identical effective stiffness tensor by exploiting a single cut section image through the original bone sample. Bone images can be extracted from CT-scan or other experimental techniques.

Following our previous works [START_REF] Hasanabadi | 3D microstructural reconstruction of heterogeneous materials from 2D cross sections: a modified phase-recovery algorithm[END_REF][START_REF] Hasanabadi | Efficient three-phase reconstruction of heterogeneous material from 2D cross-sections via phase-recovery algorithm[END_REF], several bone microstructures were generated with identical effective stiffness tensor to evaluate the effect of the geometrical distribution of bone microstructure on local energy. This enabled to study many different scenarios (different virtual samples), using the full spectrum of two-point correlation function (TPCF) for eigen microstructure, without the need of extracting real bone samples.

The characteristic function (CF) of the heterogeneous microstructure can be written [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF] as:

k n s = 1 if s 2 phase n 0 otherwise & ð1Þ
where n is a specific phase and the s index enumerates the voxel number of the sample.

This method is based on the fact that the material is not a functionally graded material (FGM) and that the statistical information of parallel cut sections are identical.

To reconstruct a three-dimensional (3D) representative volume element (RVE) based on given twodimensional cut-sections, we introduce TPCFs, being statistical functions that give some information about the volume fraction, shape and distribution of the microstructure phases and are correlated with its various features such as effective thermal and mechanical properties [START_REF] Torquato | Random heterogeneous materials: microstructure and macroscopic properties[END_REF].

In this research, we use fast Fourier transform (FFT) and a phase recovery algorithm to reconstruct the microstructure [START_REF] Fullwood | Microstructure sensitive design for performance optimization[END_REF]. The TPCF of the phases n and n 0 , C nn 0 2 , for an arbitrary vector t, can be defined as

C nn 0 2 (t) = 1 S X S s = 1 k n s k n 0 s + t ð2Þ
where S is the number of all voxels in the RVE and t is an arbitrary vector selected from the RVE. The Fourier transform (FT) of equation ( 2) leads to a more convenient form [START_REF] Fullwood | Microstructure reconstructions from 2-point statistics using phaserecovery algorithms[END_REF] for computation:

b K nn 0 m = F C nn 0 s = 1 S X S t = 1 X S s = 1 k n s k n 0 s + t e À2pi(t:m)=S ð3Þ
Using the variable exchange (t = z À s) and assuming a periodic extension of the RVE in all directions, equation ( 3) can be written as

b K nn 0 m = F C nn 0 s = 1 S X S s = 1 k n s e À2pis:m=S X S z = 1 k n 0 s e 2piz:m=S = 1 S j b K n m je Àiu n m j b K n 0 m je iu n 0 m ð4Þ
where j b K n m j and u n m are the module and phase of the FT for the CF, with

b K n m = F k n s À Á = X S s = 1 k n s e 2pis:m=S = j b K n m je iu n m ð5Þ
It should be noted that the FT, expressed as F : ð Þ, is calculated over all axes (i.e. 3D FT for (s 1 , s 2 , s 3 )). In equation [START_REF] Huiskes | Effects of mechanical forces on maintenance and adaptation of form in trabecular bone[END_REF], if n = n 0 the TPCF is also called autocorrelation function, and its FT, using equation ( 4), can be obtained as

b K nn m = 1 S j b K n m j 2 ð6Þ
Equation [START_REF] Frost | Bone ''mass'' and the ''mechanostat'': a proposal[END_REF] shows that the FT of the autocorrelation function is related only to the FT module of the CF. If the FT of the autocorrelation functions exist, it is then possible to calculate its modules. Hence, for a two-phase microstructure, if the autocorrelation functions are to be known at first, it is possible to reconstruct the microstructure if the phase of the microstructure (i.e. u n m ) is determined (see equation ( 5)). The reconstruction of an object based just on its FT modules, called phase recovery [START_REF] Fienup | Phase retrieval algorithms: a comparison[END_REF], is a common method in astronomy and microscopy.

The application of this procedure for microstructure reconstruction was first introduced by Fullwood et al. [START_REF] Fullwood | Microstructure reconstructions from 2-point statistics using phaserecovery algorithms[END_REF]. The reconstruction procedure using the phase recovery algorithm is as follow:

(1) Based on the known autocorrelation functions, the module of the CF is calculated using equation [START_REF] Frost | Bone ''mass'' and the ''mechanostat'': a proposal[END_REF]. It is called the original amplitude. From this first step a random microstructure is generated.

(2) Once the FT of the randomly generated microstructure calculated, its amplitude is replaced by the original amplitude and then the inverse FT is calculated. (3) Constraints of the generated microstructure are applied. If the difference between the FT amplitude of the obtained CF and the original amplitude is below a given threshold, the reconstruction procedure is considered completed and acceptable. Else, the obtained CF is defined inadequate for this generated microstructure and the procedure is reset until the difference between the two amplitudes falls under the given threshold.

An improvement of this method was recently proposed by Hasanabadi et al. [START_REF] Hasanabadi | 3D microstructural reconstruction of heterogeneous materials from 2D cross sections: a modified phase-recovery algorithm[END_REF][START_REF] Hasanabadi | Efficient three-phase reconstruction of heterogeneous material from 2D cross-sections via phase-recovery algorithm[END_REF] on a new approximation for the calculation of the TPCFs for a 3D RVE. We assume that the TPCF of an identical vector in a given RVE in two parallel planes are approximately equal. Hence, only two perpendicular cutsections (planes) in the RVE are necessary for the TPCFs definition and reconstruction process. This assumption is valid for statistically homogenous media [START_REF] Fullwood | Microstructure sensitive design for performance optimization[END_REF]. As shown in Figure 1, an arbitrary vector r in space can be decomposed into r z and r xy . If the TPCFs of two planes, xy and yz are known, the full set of TPCF for 3D RVE can be approximated [START_REF] Hasanabadi | Efficient three-phase reconstruction of heterogeneous material from 2D cross-sections via phase-recovery algorithm[END_REF] based on

C 11 2 r ð Þ' C 11 2 (r xy )C 11 2 (r z ) v 1 + (v 1 À C 11 2 (r xy ))(v 1 À C 11 2 (r z )) (1 À v 1 ) ð7Þ
where v 1 is the volume fraction of phase 1.

The flowchart of the reconstruction procedure is presented in Figure 2. More details of this method can be found in [START_REF] Fullwood | Microstructure reconstructions from 2-point statistics using phaserecovery algorithms[END_REF].

In this approach, TPCFs were extracted using CT images and used to reconstruct microstructures converging to the real ones. The main reasons to choose this reconstruction technique are as follows:

(i) A strong contrast expansion shows that the stiffness tensor is linked to the TPCF of the microstructure; therefore the TPCF is a very good selection of the statistical distribution descriptor. (ii) In the phase recovery algorithm, different sets of TPCF can be employed for the reconstruction procedure, therefore enabling to reconstruct a variety of microstructures with few cut-section images. An improvement of the reconstruction process is envisaged in future works using two-point cluster functions but the mathematical framework for this needs to be developed.

Finite element model

The original CT scanned images were cropped into size 351 × 319 pixels and were imported into the VCAT software (released by V-CAD Program, RIKEN, Japan). The images were then assembled into stacks of 150 images after thresholding and made a 3D RVE of the bone microstructure. Figure 3 shows the 3D RVE for the real bone specimen and the statistical one possessing similar realistic features (size, shape, and distribution) as the real one, being suitable for the calculation of its mechanical properties. Figure 3 also shows the connected network of voids inside the cortical bone. Effective properties were extracted directly from the results of the FE analysis on statistical and real bone microstructures. A direct comparison was made between the real bone sample and statistical reconstructed ones to validate the reconstruction process. In this work, the primary goal is to focus on the ability of the statistical reconstruction method to provide ''bone like'' microstructures and check whether these reconstructed microstructures show similar mechanical characteristics as the real bone. Once the statistical reconstruction method validated on a number of samples, it is then possible to create virtually an infinite number of microstructures showing different microstructures distributions and bone densities, hence enabling to study an infinite number of ''bone samples'' without the need of specimen extraction. Although the computation time is an important aspect, this is not the primary focus here as it can be addressed in future studies, but rather crucial to validate the quality of the reconstruction method in the first place.

The V-CAT software was used with its built-in mesh generator to obtain a high quality tetrahedral FE mesh, accounting for highly irregular and high aspect ratio solid and void clusters. This capability enabled to preserve a good quality of the void volume fraction during the meshing process. The mesh simplifier software has been used to reduce the number of elements in the bone part without sacrificing the volume fraction of the void. The total number of elements in the bone part was around 3.5 million tetrahedral elements for all RVEs. The volume fraction of the void phase has been calculated around 8%, which is very close to the original CT scanned data. More details about creating the 3D model, mesh generation and homogenization process for real and reconstructed microstructures are provided in [START_REF] Sheidaei | 3-D microstructure reconstruction of polymer nano-composite using FIB-SEM and statistical correlation function[END_REF].

The real internal bone microstructure either cortical or trabecular is subjected to the combination of normal and shear loads defined by the external bone geometry and load conditions (muscles, body weight, bone type). However, it is rather difficult to extract quantitative data from such conditions. Hence, in order to be able to validate the statistical reconstruction method by comparing developed stresses and elastic strain energy on real and statistical bone microstructures, we conducted six simulations (three tensile tests and three shear tests) on each RVE using the FE software ABAQUS Ò . These were used to find the modulus of elasticity and Poisson's ratios of the RVEs in the X, Y, and Z directions. Since bone is a highly heterogeneous material [START_REF] Allena | Heterogeneous directions of orthotropy in three-dimensional structures: finite element description based on diffusion equations[END_REF][START_REF] Cluzel | A general method for the determination of the local orthotropic directions of heterogeneous materials: application to bone structures using m CT images[END_REF] and this affects its overall mechanical behavior [START_REF] Altenbach | On the constitutive equations of viscoelastic micropolar plates and shells of differential type[END_REF][START_REF] Eremeyev | Material symmetry group and constitutive equations of micropolar anisotropic elastic solids[END_REF], it is necessary to compare the obtained results at the effective macroscopic scale. Based on these tests, stiffness tensors of all RVEs were obtained and influences of the microstructure distributions were quantified.

Results and discussion

Five statistical bone microstructures were reconstructed (see example in Figure 3). From these different cases, including real bone sample, homogenized macroscopic stiffness coefficients (i.e. Young modulus and Poisson's ratio) were extracted for the validation of the reconstruction process, assuming the material (bone solid phase) mechanical behavior was linear elastic with a Young modulus E = 15 GPa (average value for cortical bone) and Poisson's ratio n = 0.3. The results are presented in Table 1. The real bone sample presented values of Young modulus varying between 11.17 GPa and 12.89 GPa, and Poisson's coefficients between 0.228 and 0.301 in the different directions. It is observed that the Young moduli E x and E y are the lowest and E z is the highest which is due to the fact that cortical bone microstructure is orientated along the Haversian canal. Comparing these values with the obtained homogenized coefficients of all statistically reconstructed cases, it was observed that for all of them, the obtained Young moduli are very close to the minimum value obtained for E y of section 6 at 10.76 GPa and highest value for E z of sections 2 at 13.326 GPa. The errors developed from the microstructures variations for the Young moduli are; for the minimum Young modulus value equal to 3.1%, for the maximum Young modulus value equal to 3.38% and for the range min-max equal to 4.39%. For the Poisson's ratio, their variability checked between 0 and 0.06, which is a very small range. In addition, it was observed that the real sample mechanical properties are mostly situated in the middle of all reconstructed cases which show a good correlation with the different scenarios.

From the original bone sample, it was possible to reconstruct different bone microstructures corresponding to different case scenario (different samples/different virtual patients) with a very good accuracy. The reconstruction method allows obtaining an infinite number of virtual cases (patients) from the real bone sample, to study their mechanical behavior individually, and compare them with the real bone sample. It is also possible to extract tendencies with given scenarios of specific microstructure distributions.

Once the statistical reconstruction of different microstructure was validated, both tensile and shear mechanical loads were applied on each of them to compare their mechanical behavior, as well as the stresses and energy/energy density distribution fields developed. Each model is constituted of about 3.5 million tetrahedral FE elements in the solid phase. The applied load corresponds to an average standing person body weight (about 750 N) downscaled to the size of the sample surfaces. This leads to small deformation in the sample for the given mechanical properties; hence linear elastic analysis is adequate here. For the sake of lightness of the paper, we present only the results of the real bone sample with two other statistical microstructures built from cut section 3 and cut section 5 on the real bone sample respectively.

For the tensile test, the stresses are plotted in the corresponding loading direction (x) and (z), together with the corresponding SED developed within the microstructures. They are presented in Figures 4 and5, respectively. For the shear load scenario, stresses in the (xy) direction together with SED are presented in Figure 6. Although the effective equivalent mechanical properties are similar for all cases (see Table The interpretations of these results tend to show that:

(i) The statistical reconstruction method provides a practical method to study any virtual sample from the mechanical behavior point of view as homogenized mechanical variables (stiffness coefficients, stresses, and strain energy) remain very similar. It provides a viable method to extract local microstructures interpretations. More specifically, it shows the equivalent correlation between homogenized bone density measurements (as in magnetic resonance imaging information, for example, for the gray scale of the Hue value) and different types of microstructure distributions. (ii) The variabilities of the mechanobiological response for each reconstructed case can be evaluated as a function of the distributed SED and the different microstructure distributions. This in turn can provide long term variations in the bone density evolution and can show the influence of each of the different microstructures on the bone density kinetics.

To evaluate the potential bone density evolution as a function of the applied mechanical loads, the total developed strain energy inside each sample for the different cases (real sample and five reconstructed cases) in tensile and shear load conditions are presented in Figure 7.

Small variations of the developed strain energy are observed for each of the different load scenarios. For tensile loads, the strain energy varies from a minimum of 1.485 mJ/mm 3 up to 1.786 mJ/mm 3 . For shear loads, it varies from a minimum of 2.284 mJ/mm 3 up to 2.622 mJ/mm 3 . For all cases, the real bone sample is always around the mid-value of the statistical cases in the same way as was observed for the stiffness coefficients. The strain energy variation for all cases is within 3% to 6% compared to the real bone sample showing the reliability of the statistical reconstruction method. The developed strain energy in the z direction is slightly higher than in the x and y directions for tensile load conditions. This is in agreement with the higher stiffness of the z direction reported in Table 1. However, interestingly, shear load cases show significantly higher strain energy than simple tensile load conditions. The average energy for shear load has increased by 53% compared to tensile load conditions. This needs to be investigated in future work with a more exhaustive screening of different bone microstructure in order to extract the coupling effects between shear and tensile load conditions.

Nevertheless, as bone cells activation will depend on the developed mechanical energy within the structure, we can assume that for different types of bone microstructures, this coupling effect (between shear and tensile loads) will impact on the bone remodeling process. If shear stresses show to be predominant, it could be correlated to the fact that cells are more sensitive to shear load conditions than hydrostatic pressure [START_REF] Becquart | Human mesenchymal stem cell responses to hydrostatic pressure and shear stress[END_REF], hence favoring shear loads for optimized bone reconstruction.

In order to correlate the viability of the reconstructed microstructures to the real bone sample, probability distribution functions (PDFs) were extracted. It represents the number of elements in the sample having a specific strain energy and is plotted as an element count versus strain energy divided by the total number of elements in the model. Figure 8 shows the average PDF value of the statistically reconstructed cases with the real bone PDF distribution. Statistical variations are small and a very good agreement is observed between the reconstructed cases mechanical response and the real sample.

Since mechanical load type, load direction and microstructure distribution impact directly the bone density evolution, the SED for all the elements in each RVE were extracted. The PDFs of the SED for each load case and microstructure are presented in Figure 9.

For each load case, the PDF of real bone is plotted in red. We observe that not only the overall mechanical response of the statistically reconstructed microstructures is similar to the one of real bone but also the distributions of strain energy densities are similar. The variabilities observed in these results are based on the fact that a 3D stochastic reconstruction technique has been used for reconstruction of the bone microstructure. The proposed method includes two steps: (i) a full set of TPCFs is approximated using correlation functions being extracted from two perpendicular images; (ii) the full set of approximated TPCFs are used to reconstruct microstructures using a phase recovery algorithm. Hence, the accuracy of the reconstructing technique can directly be related to the approximation of the full set of TPCF. In our previous study [START_REF] Izadi | Effect of 2D image resolution on 3D stochastic reconstruction and developing petrophysical trend[END_REF], it was showed that for anisotropic microstructure, the error of the estimated TPCF and consequently, the error of the reconstruction procedure increased in the direction with high anisotropic distribution of microstructure. In addition, the reference image used for the reconstruction is the source of errors as it depends on the approximated TPCF in the direction perpendicular to the microstructure distribution of the reference image [START_REF] Izadi | Effect of 2D image resolution on 3D stochastic reconstruction and developing petrophysical trend[END_REF].

Nevertheless, accounting for the anisotropic complexity of the trabecular bone microstructures, the differences obtained from the reconstruction method remain small, as presented in Table 2. Comparing these results with those in Figure 9 shows that: (i) All strain energies over the whole sample for the statistical cases remain lower than 10% variations compared with the real bone sample. For simple tensile load cases, variations were between 0.2% and 4.5%, showing the influence of the anisotropic microstructure distribution over the reconstruction method [START_REF] Izadi | Effect of 2D image resolution on 3D stochastic reconstruction and developing petrophysical trend[END_REF]. For shear load cases, the differences were up to 9.4% in the yz direction showing in addition to the statistical error the sensitivity of the shear load scenario on the microstructure distribution, as presented on Figure 7. (ii) For the standard deviation variations, they are not correlated to the strain energy variations as the mean strain energy does not account for the individual microstructure distributions. However, we can still observe that for tensile cases, the z direction shows higher difference due to the stiffness difference with other directions and reconstruction method [START_REF] Izadi | Effect of 2D image resolution on 3D stochastic reconstruction and developing petrophysical trend[END_REF]. This does not appear in the shear load case as it is a tangential load. The full width at half maximum (FWHM) deviations of the PDF for tensile load remain within the 10% variability. For the z direction, it is almost identical. It appears from the results in Table 1, Figures 7 to 9, and Table 2 that the statistical reconstruction method enables to study virtually an unlimited number of samples (bone like microstructures) as the calculated effective mechanical properties are very similar to the one of real bone. Some improvements remain to be done as the distribution of the strain energy shows variabilities over the statistical microstructures (up to 10%) due to the abilities of the reconstruction method to best adapt to the real bone microstructure. These variations need to be evaluated with their respective impact on the bone density evolution over long periods of time in future works.

i) (ii) (iii) (iv) (i) (ii) (iii) (iv) (i) (ii) (iii) (iv) (i) (ii) (iii) (iv) (
Figure 10. Strain energy density (mJ/mm 3 ) on elements surrounding the void space for mechanical load applied in the x, y, z, and xy directions: (i) real bone, (ii) cut section 3, (iii) cut section 5, and (iv) probability distribution function for the three cases.

As bone remodeling occurs on internal surfaces of bone where osteoblasts and osteoclasts are present, it is necessary to look at the developed strain energy densities on these surfaces. These are presented on Figure 10 for three cases : (i) real bone, (ii) cut section 3, and (iii) cut section 5 for the mechanical load conditions of simple tension in x, y, and z directions and shear xy. In addition, the PDF functions vs SED are plotted for the elements on the border. The tail of the PDF plots beyond a given energy is not shown on the figures as their PDF values are very small.

The presented results highlight the geometrical localization where bone remodeling is more likely to occur as osteoblasts and osteoclasts (respectively responsible of construction and resorption) are located on the bone surface. Similar to previous results, the maximum distributed SED is for shear load scenario with an average value around 1 mJ/mm 3 whereas for tensile loads, the average values are around 0.2 mJ/ mm 3 for x and y tensile loads and 0.7 mJ/mm 3 for z tensile load. It appears again that tensile loads provide similar results compared to the stiffness coefficients proportions. We also observe that the real bone sample shows mostly higher SED for tensile loads but lower for shear load. The statistical samples are very close to each other, which could be correlated with the reconstruction method, as explained above [START_REF] Izadi | Effect of 2D image resolution on 3D stochastic reconstruction and developing petrophysical trend[END_REF].

Nevertheless, as far as the mechanical properties are concerned, a very good correlation is observed between the statistical cases and real bone which therefore provides good confidence in the statistical distribution of the developed strain energy. The impact of the observed differences needs to be quantified with regard to the biological distribution of cells on the developed surface microstructures and cell activation thresholds for bone remodeling. This will be done in future works.

Conclusion

A statistical reconstruction method was utilized to reconstruct cortical bone-like microstructures. The method is based on TPCFs using two perpendicular cutting planes within a real bone sample. Unlimited statistical bone microstructure distributions can be obtained with this method enabling the study of the influence of bone microstructure distribution without the need of extracting many patient bone samples.

All statistically reconstructed samples had identical effective bone density compared to the real bone sample and showed similar mechanical properties. The local stress, strain, and strain energy were calculated for all samples under simple mechanical load conditions. Local variations were observed as a function of the different reconstructed microstructures, and linked to effective values. Good correlation was observed for all cases.

The statistical reconstruction method enables us to create virtual bones, extract the variations of the local strain energy distribution as a function of the bone microstructure, and compare with macroscopic one. From this, bone remodeling could be incurred as a function of the local cell activation processes on the bone inner free surface. This local bone remodeling could therefore be linked to the effective bone density. The knowledge of the anisotropic stiffness tensor could then be used as a complementary indicator of the long term evolution of bone density. This needs to be validated in future work.

Figure 1 .

 1 Figure 1. Decomposition of the vector r to r z and r xy .

Figure 2 .

 2 Figure 2. Flowchart of the reconstruction procedure.

Figure 3 .

 3 Figure 3. Finite element mesh for cortical bone microstructures (left: bone phase, right: void phase): (a) real specimen; (b) statistically reconstructed specimen from cut section 2. The approximate size of samples (real and reconstructed ones) is 1930 x 1750 x 600 mm 3 .

1 )

 1 , the stresses and SED, presented in Figures4 to 6, show different distributions for both tensile and shear applied loads. Despite this anisotropic distribution, the amount of developed energy within the samples remains similar (see below).

Figure 4 .

 4 Figure 4. Developed stresses and strain energy density in the real bone and statistical samples (section 3 and 5) under tensile load applied in the x-direction.

Figure 5 .

 5 Figure 5. Developed stresses and strain energy density in the real bone and statistical samples (section 3 and 5) under tensile load applied in the z-direction.

Figure 6 .Figure 7 .

 67 Figure 6. Developed stresses and strain energy density in the real bone and statistical samples (section 3 and 5) under shear load applied in the xy-direction.

Figure 8 .

 8 Figure 8. Probability density function for the YY (a) and ZZ (b) tensile cases, and XY (c) shear case as a function of the developed local strain energy (mJ/mm 3 ).

Figure 9 .

 9 Figure 9. Normalized probability distribution function of strain energy density for the different microstructures with different mechanical load directions.

Table 1 .

 1 Modulus of elasticity and Poisson's ratio in different directions for the real and statistical bone microstructures.

		E x (GPa)	E y (GPa)	E z (GPa)	n xy	n xz	n yz	n yx	n zx	n zy
	Real bone	11.723	11.172	12.894	0.228	0.263	0.262	0.275	0.289	0.301
	Cut section 2	12.506	12.064	13.326	0.291	0.272	0.274	0.281	0.29	0.301
	Cut section 3	11.456	11.310	12.686	0.281	0.260	0.264	0.277	0.288	0.298
	Cut section 4	11.456	11.310	12.713	0.281	0.260	0.264	0.277	0.289	0.297
	Cut section 5	11.556	11.229	12.838	0.288	0.259	0.263	0.279	0.288	0.297
	Cut section 6	11.185	10.765	12.370	0.282	0.259	0.258	0.271	0.287	0.298

Table 2 .

 2 Differences between real bone sample and averaged of statistical cases for: (i) strain energy density variation in different directions in tensile and shear load condition (X, Y and Z indicate tensile test and XY, XZ and YZ indicate shear tests), and (ii) standard deviation difference at full width and half maximum (FWHM) as a representative of the variation of the energy density of each element compared to the mean value.

		Tension X	Tension Y	Tension Z	Shear XY	Shear XZ	Shear YZ
	Mean strain energy density variation (mJ/mm 3 )	0.0186(1.8%)	0.0449(4.5%)	0.0021(0.2%)	0.0699(7%)	0.0138(1.38%)	0.0942(9.4%)
	Standard deviation	11%	10%	42%	6%	11%	7%
	difference at FWHM

Funding

This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.