Kafil Hajlaoui
email: kafil.hajlaoui@inist.fr

Pascal Cuxac
email: pascal.cuxac@inist.fr

Jean-Charles Lamirel
email: lamirel@loria.fr

Claire François
email: claire.francois@inist.fr

Jean Charles Lamirel

Enhancing Patent Expertise through Automatic Matching with Scientific Papers

Keywords: Supervised classification, Technological and scientific survey, Patents, KNN, Association rules

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Text categorization is a machine learning task which aims at automatically assigning predefined category labels to new upcoming free text documents with related characteristics [START_REF] Cohen | A survey of current work in biomedical text mining[END_REF]. Because of its numerous applications, text categorization has been one of the most studied branches within the field of machine learning [START_REF] Hillard | An active learning framework for classifying political text[END_REF]. Consequently, a variety of classification algorithms were developed and evaluated in applications such as mail filtering [START_REF] Cormack | Online supervised spam filter evaluation[END_REF], opinion and feelings analysis [START_REF] Pang | Opinion mining and sentiment analysis[END_REF], news [5] [6] or blogs [START_REF] Durant | Predicting the Political Sentiment of Web Log Posts Using Supervised Machine Learning Techniques Coupled with Feature Selection[END_REF] classification. Among the most often used learning methods exploited in that context, we may mention artificial neural networks [START_REF] Wiener | A Neural Network Approach to Topic Spotting[END_REF] [9], K-nearest-neighbors (KNN) [START_REF] Yang | An example-based mapping method for text categorization and retrieval[END_REF], decision trees [START_REF] Lewis | Comparison of two learning algorithms for text categorization[END_REF] [12] [START_REF] Apte | Text mining with decision rules and decision trees[END_REF], Bayesian networks [START_REF] Lewis | An Evaluation of Phrasal and Clustered Representations on a Text Categorization Task[END_REF] [15], support vector machines (SVM) [START_REF] Joachims | Text categorization with support vector machines: Learning with many relevant features[END_REF], and more recently, boosting based methods [START_REF] Schapire | Boosting and Rocchio applied to text filtering[END_REF] [START_REF] Iyer | Boosting for document routing[END_REF]. Although many methods developed for automatic text categorization have achieved significant accuracy when applied to simple text structure (for example emails, summaries, etc.), there are still many remaining challenges concerning classification of complex documents, especially when classification relies on imbalanced learning data.

A broad range of studies address the problem of Medline2 database categorization. Most of these works focus on the importance of data preprocessing and data representation steps in the context of the text categorization task. In [START_REF] Lan | Text representations for text categorization: a case study in biomedical domain[END_REF], the authors show that, in the case of a text representation based on the "bag of words" model, the weighting of the extracted terms significantly increases the performance of the classifiers. In order to classify Medline papers into predefined topics, Suomela and Andrade [START_REF] Suomela | Ranking the whole MEDLINE database according to a large training set using text indexing[END_REF] restrict the extracted descriptors to predefined lexical classes (nouns, adjectives, verbs) and apply a word frequency scheme. Using specific Medline topics, the authors obtain a classification F-score of 65%. The same approach is further used by the Medline Ranker web-service [START_REF] Fontaine | MedlineRanker: flexible ranking of biomedical literature[END_REF] whose role is to extract a relevant list of Medline references starting from a set of keywords defined by the user. The study of Yin et al. [START_REF] Yin | Document classification for mining host pathogen protein-protein interactions[END_REF] focuses both on the identification and on the extraction of protein interactions from Medline articles. For that purpose, documents are preprocessed using bigrams, in a first step, and SVM method is applied, in a second step. The authors obtain a performance of 50% true positives, and a recall rate of 51%. Recently, the Biocreactive III challenge proposed to classify Medline articles belonging to the biomedical field [START_REF] Krallinger | Results of the BioCreative III (Interaction) Article Classification Task[END_REF]. The best performance related to Medline data was obtained on this collection, with an accuracy of 89.2% and an F-score of 61.3%.

Up to now, patents evaluation is a manual operation that involves groups of experts with in-deep expertise of the related field. This evaluation is mainly based on references to relevant scientific papers (articles, theses, books ...) associated to the patents. The automatic classification of publications in patent classes can thus represent a valuable help for the experts. However, this task is not a traditional classification task because the classification structure (i.e. the patent classification scheme) does not directly fit with the data to be classified (i.e. the scientific papers). To cope with that problem, two alternative approaches can be used. A first approach consists in creating a gateway between the publications classification plan and the patents classification plan. Nonetheless, this approach is difficult to implement because it involves the intensive use of complex tree comparison techniques (here, the classification plans), and consequently, an intensive use of complementary human expertise. A second approach consists in developing a classification system that directly uses the patents classification plan. Such an approach is founded on the assumption that scientific papers that are cited in a patent are strongly related to the patent field, and consequently to the classification code of the latter. In this context, the training dataset would consist of the whole set of scientific citations extracted from the patents of the considered field. However, one potential barrier of this approach is that the learning classes might not necessarily have a homogeneous quantity of patents and thus might not provide an homogeneous amount of learning data (i.e. cited papers) leading thus to face with an imbalance classification problem. Moreover, in a focused domain, patent classes might have high citation or cited topics overlaps leading to additional class similarity problem.

In the following sections, we describe a complete experiment of automatic classification of scientific papers based on an initial patents dataset. The experimental dataset is related to the field of pharmacology and the bibliographic references cited in the patents are extracted from the Medline collection. In the first section, we describe the dataset construction process and we illustrate the resulting phenomena of imbalance learning examples and class similarity. By applying usual categorization methods, we then illustrate, in the second section, the influence of the term extraction strategy on the classification results. Two approaches are particularly discussed. The first one is based on the direct use of Medline indexing keywords associated to the bibliographic records. The second is based the construction of an index from the titles and abstracts of the records by the use of NLP tools. This section highlights that the best performance are obtained with the K-Nearest Neighbor algorithm (KNN) in our context. To cope with the class imbalance and similarity problems, we present in the third section a modified KNN algorithm named KNNBA-2T which is based on the exploitation of association rules between data descriptors and patent class labels. We show that this algorithm provides better classification accuracy than the original KNN algorithm in our context. In section 4, we perform a complementary test of the KNNBA-2T algorithm in combination with resampling techniques. In this test, we exploit our former dataset as well as 6 other UCI datasets and compare the results of KNNBA-2T with resampling with a broader range of other usual algorithms. The final section draws our conclusion and perspectives.

Building and indexation of the corpora

Our main experimental resource is issued from the QUAREO project. It is a collection of patents related to the pharmacology domain and including bibliographic references. This resource consists of 6387 patents in XML format, grouped into 15 subclasses of the A61K class (medical preparation). As shown in Figure 1, we begin by extracting the citations from the patents. From 6387 patents, we extracted 25887 references such as databases, books, encyclopedia and scientific articles. In a second step, we query the Medline database with the extracted citations related to the scientific articles. In such a way, we obtained 7501 articles. This represents a recall of 90% for this type of references. Each article is then labeled by the class code of the citing patent and the set of labeled articles represents our final training dataset.

Data Representation

As we have mentioned it before, for text classification, the choice of a document model is a crucial step. A common approach is to use a document model called "bag of words", in which the only exploited information is the presence and/or the frequency of terms. In our case, we translate the bag of words model into a vectorial representation, as it has been proposed by Salton [START_REF] Salton | Automatic processing of foreign language documents[END_REF]. Following this approach, each article of the dataset is represented as a vector in an N-dimensional space, where N is the total number of terms (features) issued from the articles collection. The whole text collection is then represented as a (N +1) J matrix, where J is the number of articles in the collection. Each line j of this matrix is an N-dimensional bag of word vector for the article j, plus its class label.

If a feature i does not occur in the article j, then the relevant matrix element aij is zero, otherwise it is assigned a positive value or weight. The way to calculate this weight depends on the scheme used for feature representation. The weight is 0 or 1 for the binary scheme. On its own side, the standard frequency weighting scheme is based on terms document occurrences. However, with such scheme, too much importance could be given to descriptors that appear frequently in many documents and which are, consequently, unrepresentative for each single document. Another weighting scheme, called the TF.IDF (Term Frequency Inverse Document Frequency) is thus often used in literature [START_REF] Vincarelli | Indexation de documents manuscrits[END_REF] [27] [START_REF] Spark-Jones | A statistical interpretation of term specificity and its application in retrieval[END_REF]. This scheme evaluates the importance of a term according to its frequency in the document (TF = Term Frequency) weighted by its frequency in the corpus (IDF = Inverse Document Frequency).

. ( , ) ( , ) ( )

k j k j k Tf Idf t D TF t D Idf t = × where ( , ) k j TF t D is the number of occurrences of k t in j D , and, ) ( log ) ( k k t DF S t Idf =
where S is the documents number in the corpus and ( )

k DF t is the number of documents containing k t .
At the following stage, we built features according to two different approaches, the first one relying on keywords found in documents, and the second one relying on the lemmas extracted from the document abstracts by the use of an NLP tool. The objective of this last approach is to improve the representation of documents' content. To do this, we use the TreeTagger tool [START_REF] Schmid | Probabilistic part-of-speech tagging using decision trees[END_REF] developed by the Institute for Computational Linguistics of the University of Stuttgart. This tool is both a lemmatizer and a tagger. A lemmatizer associates a lemma, or a syntactic root, to each word in the text and a tagger automatically annotates text with morpho-syntactic information. In our case, the document are firstly lemmatized and the tagging process is performed on lemmatized items (in the case when a word is unknown to the lemmatizer, its original form is conserved). The punctuation signs and the numbers identified by the tagger are deleted. A sample output of the TreeTagger program is given in figure 3. The feature selection according to grammatical categories allows identifying salient features for the documents classification according to document types or opinions. Moreover, this approach permits to consistently reduce the description space. We thus choose to apply it in our experiment.

The

Classification

To evaluate the relevance of different indexation and weighting methods, we chose to use three different classifiers: a K-nearest-neighbors (KNN) classifier relying on Euclidean distance, a (SVM) classifier and a probabilistic classifier (Naive Bayes). These three supervised machine learning algorithms are known to provide the best results for text classification [START_REF] Sebastiani | A tutorial on automated text categorisation[END_REF] [START_REF] Yang | A reexamination of text categorization methods[END_REF]. In our case, we have exploited them in the Weka environment 4 .

In all tests, we have applied different weighting techniques, according to type of extracted descriptors. For the lemmas we mainly use the standard frequency and the TF.IDF techniques. Applying the TF technique on keywords would be meaningless, because the indexing on documents is not redundant. Therefore we use solely the Boolean or IDF technique in this case.

For the features based on lemmas, we have performed several experiments by switching the selected grammatical categories (A: Adjective, N: Noun, NA: Noun + Adjective, NV: Noun + Verb, VA: Verb + Adjective, NVA: Noun + Verb + Adjective).

The classification results are expressed in terms of precision and recall. A precision of 100% means that all articles are classified in the correct category. The recall is the percentage of correct answers that are given. These measures are calculated after applying a 10-fold cross-validation process (90% of the corpus is used for learning and 10% for testing).

Tables 1 and2 show the obtained results in terms of precision (P) and recall (R). They illustrate that the best results in our corpus are obtained with the KNN method combined with an indexation based on the lemmas involving the three grammatical categories (Nouns, Verbs, Adjectives) and using TF-IDF weighting scheme. The obtained measures are 61% for precision and 55% for recall. However, these results can still be considered as far from satisfying ones. Such results can be explained by the imbalanced data distribution between the classes (see figure 2), but also by the fact that the classes are very similar one to another. To highlight that problem, we computed class/class similarity using cosine correlation and drew the resulting class/class similarity distribution (figure 4). This distribution clearly shows that it might be difficult for any classification model to precisely detect the right class: more than 70% of classes' couples have a similarity between 0.5 and 0.9.

Class pairs count

Cosine similarity

Fig. 4. Class to class similarity distribution

We therefore propose, in the next section, an improvement of the approach, based on the best method, namely the KNN method. The goal of this improvement is to take into account the specific characteristics of the corpus: the important imbalance between classes and the high similarity between them.

4

The KNNBA-2T method

Our improvement of the KNN algorithm is based on the exploitation of association rules. We firstly present a general definition of association rules. We then present a new approach for calculating the weight of class attributes or features, by using a special type of association rules. We finally present a new algorithm, called KNNBA-2T, inspired by the method previously developed by Mordian et al. [START_REF] Mordian | k-Nearest Neighbor Based Association Algorithm[END_REF].

Association rules

The association rules extraction approach is a method for discovering relevant relationships between two or more variables. This method is based on local laws and requires no user intervention (it lets the system self-organizing). It allows identifying, from a set of transactions, a set of rules that express a possibility of association between different items (words, attributes, concepts …). A transaction is a series of items expressed in a given order. In addition, transactions can be of different lengths. The relevance of a rule of such extracted association is measured by its index of support and its index of confidence.

For an association rule: X Y → the indices of support and confidence are defined by the following two equations:

), (

Y X P Support U = ) ( X Y P Confidence = where ) ( Y X P U
is the probability that a transaction contains both X and Y, and is the conditional probability of Y knowing that it is X.

The first efficient method for extracting such rules was introduced by Agrawal for the analysis of the market basket through the Apriori algorithm [START_REF] Agrawal | Fast algorithms for mining association rules in large data bases[END_REF]. The operation of this algorithm can be decomposed into two phases:

1) Searches for all the "patterns" or frequent itemsets that appear in the database with a frequency greater than or equal to a threshold defined by the user, called minsup. 2) Generation, from these common patterns, of all the association rules with a measure of confidence greater than or equal to a threshold defined by the user, named minconf.

The KNNBA-2T algorithm

KNNBA is an improvement of the KNN algorithm. The objective is to assign weights to each attribute by using association rules. For that purpose, we used the association rules that help to identify the most representative terms of a given class. Each transaction consists of all of the extracted terms (attributes) and the label of the class. After the generation of the rules, we are keeping the rules of the type:

Class Attribute → and Class Attribute Attribute → 2 1 ,
The rules composed of three attributes are rare and thus not determinative.

The principle behind this approach is that if two attributes (here Attribute1 and At-tribute2) are associated together to a class, the relevance (i.e. the information power) of each of the two attributes deducted from their association must be considered as more important than the one of each single attribute.

The first version of our algorithm (KNNBA-1T) is similar to the Mordian et al. algorithm [START_REF] Mordian | k-Nearest Neighbor Based Association Algorithm[END_REF]. It takes into account the rules composed of a single attribute (term). In the second version of our algorithm (KNNBA-2T), we first compute two attributes rules and we further apply the former principle by deriving single attribute rules from two attributes rules.

After the rule extraction step, a weight can be associated to each attribute i (i.e. feature) [START_REF] Mordian | k-Nearest Neighbor Based Association Algorithm[END_REF]. It is computed as:

[ ] 1 [ ] ( ) 1 sup i W i G = -
where G_sup represents the greatest support of the attribute i.

As compared to KNN, the new formula for calculating the distance used in the KNNBA2T takes the weight of the attributes into account and thus becomes:

∑ = - × = n i bi ai x x i W b a D 1 2 ) ( ] [ ) , (
where a and b are two documents, and x ai and x bi represent the term i of each document vector.

Fig. 5. General operating process of the KNNBA approach

The general process of KNNBA-2T approach is summarized in figure 5 and is composed of three stages: Stage 1: this phase consists of two steps. The first step is the construction of transactions representing entries to generate the association rules. Each document is transformed into a transaction, consisting of all the representative document descriptors associated with the label of the class. The second step is the generation of the rules of association through using an Apriori algorithm [START_REF] Cohen | A survey of current work in biomedical text mining[END_REF].

Stage 2: in this phase, we seek to generate an attribute weight vector from the description of the documents. For that purpose, a group of 15 rules (15 corresponding to the number of classes) is built for each attribute and the most relevant rule (the highest support, the highest confidence) is used to compute the attribute weight.

Stage 3: this phase consists in applying the KNN algorithm with the added extension. To predict the class of a new document by the inter-document similarity calculation, we take into account the weight vector generated in the previous stage.

From an overall perspective, our approach extends the K-nearest-neighbors method in two ways:

1) First, a new weighting scheme of descriptors is introduced, according to their informational weight in relation with their distribution in all the classes.

2) Second, the vote of the closest neighbors is based on a vote function extended by the weight vector W. This extension uses the strength for each term to activate the classes.

Our proposed extension is thus founded on the general idea that the observations for training data, which are particularly close to the new observation (y, x) from test data, must have a higher weight in the decision that the neighbors that are farther from the pair (y, x). This strategy differs from for the standard KNN method in which only the K-nearest-neighbors influence the prediction, but the influence is identical for each of the neighbors, irrespective of their degree of similarity with (y, x). In order to achieve our new goal, the distances on which the search for neighbors is based on are thus transformed in a first step according to the strength (i.e. power) of the term to activate a class. 3 illustrates the precision and recall results obtained after application of the three KNN algorithms (KNN, KNNBA-1T, KNNBA-2T) on our reference dataset with the use of NVA lemmas and TF-IDF weighting. The best results are obtained with the KNNBA-2T algorithm, when compared to the KNN and the KNNBA-1T algorithms. We find that the percentage of correctly classified documents rises from 61% to 65% with KNNBA-1T and to 67% with KNNBA 2T. Our adapted methods thus significantly improve the classification performance on our dataset. Fig. 6. Correction of the class imbalance and class similarity using the KNNBA-2T approach Our new KNNBA-2T approach provides correction capabilities both for terms distribution within classes (i.e. class imbalance) and for class similarity. The correction of class similarity is illustrated at figure 6. However, we also remarked that the smoothing of terms distribution is not effective on the largest class (A61K31) which is always a majority class. As a result, our approach still tends to ignore small classes while concentrating on classifying the larger ones accurately.

Even if the class similarity problem remains difficult to solve in our context, because of the overlapping nature of the exploited patent classification, one complementary approach can be used to better solve the class imbalance problem. Hence, resampling methods are very commonly employed for dealing with such problem. Their advantage over other methods is that they are external and thus easily portable and very simple to implement. Resampling is usually conducted using the two following strategies: oversampling consists of copying existing training examples at random and adding them to the training set until a full balance is reached. Undersampling consisted of randomly removing existing examples until a full balance is reached [START_REF] Tomek | Two modifications of CNN[END_REF] [START_REF] Kubat | Addressing the curse of imbalanced training sets: onesided selection[END_REF].

In the next section we thus provide an extended test of the KNNBA-2T algorithm by combining it with a resampling technique. In this test, we exploit our former dataset as well as 6 other UCI datasets and compare the results of KNNBA-2T with resampling with a broader range of other usual algorithms.

Extended experimental results

In this new experiment we have made use of a combination of the KNNBA-2T method with a resampling technique. The exploited resampling technique is the Weka resample algorithm which is an undersampling algorithm suitable for decreasing the influence of very large classes. We also extend the range of comparison by using a broader range of classification techniques (including neural network (ANN) and J48 algorithms). Apart of our former dataset, we also exploit 6 complementary reference datasets issued from the UCI machine learning database collection. The overall characteristics of the experimental datasets are presented in table 4. All algorithms are executed in the following similar conditions:

1. The Weka resample (undersampling) algorithm is applied.

2. The ten-fold cross validation method is used.

3. The Weka default parameters are used for all of the algorithms.

4. The number of neighbors is set to 10 (k=10) for the KNN-based algorithms. In all experiments, the accuracy of each algorithm is based on the percentage of correctly classified documents. The complete results are presented in the table 5.

Our new experiments highlight that the resampling method significantly improves the performance of the KNNBA-2T method on our reference dataset of scientific papers (+10 points of precision). Even if it is not presented in the table, similar improvement can been observed for the other methods and for the other datasets. Table 5 also highlights that the overall results of our KNNBA-2T algorithm are above average (most of the time the best) on the other UCI datasets. However, the most important advantage as regards to the other methods is observed for large test collections in which classes include a large number of attributes. In this case our method clearly allows reducing the class representation space by selecting the relevant attributes.

Discussion and Conclusion

The classification of scientific papers in a patents' classification plan is a real challenge as such classification plan is very detailed and not very suitable with respect to the scientific documents. In this paper we presented a new method for supervised classification derived from the KNN method. This method, which we named KNNBA-2T, operates a classes' descriptor term weighting based on association rules induced by these terms. We applied it on a dataset of bibliographic articles from the Medline database in order to classify them within a classification plan of patents belonging to the field of pharmacology. This new method offers very interesting performance for our study as compared to existing methods, especially when it is combined with resampling techniques. Nevertheless, the resulting class imbalance and the similarity of the class description remains a major problem still hampering the performance improvement of automatic classification of articles within the international patents' plan. Therefore, we undertook new experiments in order to combine our method with vocabulary extension techniques based on domain ontologies. In our context, such ontology as Mesh which is associated to Medline resource represents a good candidate.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Building steps for the training dataset

Table 1 .

 1 Classification results related to indexing by keywords

		KNN			NB				SVM
	Boolean	IDF		Boolean	IDF		Boolean	IDF
	P	R	P	R	P	R	P	R	P	R	P	R
	0.39	0.39	0.39	0.43	0.4	0.47	0.43	0.44	0.4	0.45	0.4	0.45

Table 2 .

 2 Classification results related to indexing by lemmas

			KNN				NB				SVM
		Frequency	TF.IDF	Frequency	TF.IDF	Frequency	TF.IDF
	Type	P	R	P	R	P		R	P	R	P		R	P	R
	A	0.42 0.36 0.42 0.36 0.38 0.2	0.37 0.18 0.45 0.46 0.45 0.46
	N	0.5	0.41 0.52 0.4	0.43 0.31 0.44 0.28 0.54 0.55 0.54 0.55
	NA	0.55 0.4	0.57 0.39 0.45 0.36 0.46 0.36 0.55 0.55 0.55 0.55
	NV	0.49 0.38 0.52 0.38 0.44 0.35 0.44 0.31 0.53 0.54 0.53 0.54
	NVA 0.6	0.54 0.61 0.55 0.44 0.34 0.45 0.34 0.54 0.55 0.55 0.55
			40
			35
			30
			25
			20
			15
			10
			5
			0
				0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1

Table 3 .

 3 Comparison of classification results with KNN and KNNBA algorithms

	KNN	KNNBA-1T	KNNBA-2T
	P	R	P	R	P	R
	0.61	0.55	0.65	0.65	0.67	0.67
	Table

Table 4 .

 4 Description of the datasets used in the experiments

		Dataset	Size	Nb. of attributes	Nb. of classes
	1 NVA+Resample	7501	2463	15
	2 Breast-cancer-w	699	11	2
	3	Car	1728	6	4
	4	Ecoli	336	6	8
	5	Glass	214	10	6
	6	Nursery	266	8	5
	7	Zoo	101	18	7

Table 5 .

 5 Comparison of accuracy of KNNBA-2T with other algorithms

	Dataset	KNNBA-2T KNN	NB	NN	J48	SVM
	NVA+Resample	77.78	70.2	70.1	71.56	68.27	71.56
	Breast-cancer-w	96.89	96.42	95.99	95.27	94.56	96.99
	Car	95.21	93.51	85.53	99.53	92.36	93.75
	Ecoli	88.31	86.01	85.41	86.01	84.22	84.22
	Glass	68.89	66.35	48.59	67.75	66.82	56.07
	Nursery	98.58	97.58	95.08	99.83	98.06	96.93
	Zoo	98.41	88.11	95.04	96.03	92.07	96.03

http://www.quaero.org

http://www.ncbi.nlm.nih.gov/pubmed/

http://www.cs.waikato.ac.nz/ml/weka/index.html

http://www.quaero.org

http://www.oseo.fr/

Acknowledgment: This work was done under the program QUAERO 5 supported by OSEO 6 national agency of research development.


Ivana Roche
email: ivana.roche@inist.fr

Nathalie Vedovotto
email: nathalie.vedovotto@inist.fr

Claire François
email: claire.francois@inist.fr

Dominique Besagni
email: dominique.besagni@inist.fr

Marianne Hörlesberger
email: marianne.hoerlesberger@ait.ac.at

Dirk Holste
email: dirk.holste@ait.ac.at

Edgar L Schiebel
email: edgar.schiebel@ait.ac.at

Pascal Cuxac
email: pascal.cuxac@inist.fr

Assessment of the applied orientation of a researcher's production: An informetric approach based on a content analysis 1

Assessment of the applied orientation of a researcher's production: An informetric approach based on a content analysis

Extended abstract

One way of making the distinction between pure and applied research was introduced by Donald Stokes [START_REF] Stokes | Pasteur's Quadrant -Basic Science and Technological Innovation[END_REF] who defined a two dimensions chart, "the Pasteur's Quadrant". It is a label given to a class of scientific research developments that both seek fundamental understanding of scientific problems, and, at the same time, seek to be eventually beneficial to society. The works of Louis Pasteur, a French chemist and physicist, pioneer of the microbiology, are thought to exemplify this type of method, which bridges the gap between "basic" and "applied" research. The Pasteur's Quadrant characterizes three distinct classes of research:

-pure basic research, illustrated by the work of Niels Bohr, early 20th century atomic Danish physicist; -pure applied research, exemplified by the work of Thomas Edison, North-American inventor and businessman; -use-inspired basic research, described as "Pasteur's Quadrant".

A classic solution to determine the applied orientation of the works of a researcher is to search for the possible patents in the development of which he(she) took part (e.g. [START_REF] Glänzel | Patents cited in the scientific literature: An exploratory study of "reverse" citation relations[END_REF][START_REF] Moed | Handbook of quantitative science and technology research: The use of publication and patent statistics in studies of S&T systems[END_REF][START_REF] Glänzel | Publication activity, citation impact and bi-directional links between publications and patents in biotechnology[END_REF]. Usually, the submission of a patent is the outcome of an initiative of technology transfer and can be indeed considered as the practical fulfilment of research works presenting determinedly applied characteristics. Another possible approach consists on directly examining the researcher's works published in the S&T (scientific and technological) literature and characterizing their content in applied or fundamental.

In the context of a previous European project (Holste, 2012), we modelled this two criteria by means of indicators calculated as:

-the enumeration of the patents in which the researcher has contributed. This indicator gets integer values included in the interval [0,  [. Unfortunately, the number of patents is often very low, which involves a lack of accuracy of the related indicator.

1 The authors wish to acknowledge the support of the IDEAS specific programme of the EU's 7th Framework

Programme for Research and Technological Development in which context this work is partially accomplishing (DBF project -reference n° 240765 -within the Coordination and Support Actions).

Page 2 of 3 -the ratio of the researcher publications appearing in journals which content is categorized as applied. This indicator gets real values between 0 and 1. For both indicators we guessed that the higher the value, the more asserted is the applicability of the researcher works.

Even if such an approach could sound pragmatic, it conceals a few weaknesses in particular because of the journals categorization step. Indeed, by using a binary categorization "applied or fundamental", it seems a priori easy to automatically transpose the journal's category to all the articles that are published in it. A first difficulty appears when we have to determine the criteria defining this binary categorization of the S&T journals. If works have been done supplying hierarchical tree classifications, more or less detailed, of the scientific domains, the problem persists when we wish to determine which ones can be determinedly considered applied or fundamental.

In addition, the category of a journal can be not unique but variable, according to the scientific domain of each researcher that publishes his(her) work in it. Let us consider, for instance, the Biology domain as a priori fundamental. All the journals classed in this domain get then the category "fundamental", as well as all the articles published in them. Yet this does not always correspond to reality. Indeed, if this remains true for the biologists' publications in these journals, that of an IT specialist who would bring a development software to Biology should receive the "applied" category.

And so in this work, instead of examining directly the researcher's production, we propose to assess its applied orientation by analyzing the S&T literature citing his(her) publications. The citation expresses its exploitation, in different ways and at different degrees of importance, and we consider that this corpus of bibliographic records referring to the researcher's production represents a real and pragmatic information source about the utilization of his(her) former works by the scientific community in new developments. The content analysis approach applied to this corpus gives us the means to appreciate the degree of application of the researcher's scientific production. This way, we can detect potentially applicable works whose results could be integrated in more applied issues. So, in a first step, we extract from a bibliographic database the list of the researcher's publications in S&T literature. This list enables us to determine the set of publications citing at least one of the extracted publications. This corpus can be considered as an image of the scientific landscape of citing papers that are based on the past work of this researcher.

Then a data mining step, based on NLP (natural language processing) techniques, is used to obtain an assisted indexing of the records by combining them with keywords. A clustering is then applied to this enriched corpus. The software tool uses a non-supervised and non-hierarchical clustering algorithm, the axial K-means, inspired by Kohonen's self-organizing maps formalism. This method considers the keywords as indicators of the content of bibliographic records, which in their turn are considered as indicators of the research themes. This first step is followed by a principal component analysis leading to a 2D-mapping of the clusters. Thematic networks emerge from the relations between clusters and, according to a geographical metaphor, build a map of the corpus research landscape.

At this stage, a scientific expert performs an analysis of the clusters in terms of their content and of their position and relations in the map. For this analysis, he(she) must adopt a particular point of view, in order to evaluate whether the content of a corpus is essentially applied or essentially fundamental: he(she) looks therefore at the content of each cluster by considering the bibliographic records title and keywords, in order to evaluate how much the content of the clusters can be considered as applied. We call this the level of "applicateness" of the cluster. Finally, we weight the "applicateness" value of each cluster with the help of a parameter coming from the clustering results, in order to take into account the global context of the corpus research landscape represented by the map of thematic clusters and, in particular, by analyzing the links between them. The combination of the weighted "applicateness" of each cluster gives us the level of the applied orientation of the researcher's production.

In this paper, we will firstly detail the approach used to evaluate the applied orientation of a researcher's production. Secondly, we will describe the data processing, and finally we will, on one hand, present and discuss the obtained results for a few case studies and, on the other hand, look towards the future by proposing possible improvements of our methodology.


Régis Clouard

Abderrahim Elmoataz

Marinette Christine Porquet


M Saint-Paul
email: michel.saint-paul@neel.cnrs.fr

C Guttin

A Abbassi

Zhao-Sheng Wang

Huiqian Luo

Xingye Lu

Cong Ren

Hai-Hu Wen

K Hasselbach

Surface impedance of BaFe 2-x Ni x As 2 crystals

Keywords: A. Pnictides, D. Superconductivity, E. Surface Impedance,

Measurements of the real  1 and imaginary   part of the conductivity were performed in optimally doped BaFe 1.9 Ni 0.1 As 2 and overdoped BaFe 1.88 Ni 0.12 As 2 crystals in the frequency range 20MHz -1.5 GHz using a single coil technique. The temperature dependence of the London penetration depth follows a T 2 law. The conductivity  1 increases with decreasing temperature below Tc in agreement with the results obtained for the optimally Co doped BaFe 2-x Co x As 2 crystals. The increase of  1 in the superconducting state is attributed to a rapidly decrease of the quasiparticle scattering rate.

Introduction

A number of experimental works have been published recently concerning the surface impedance of iron based superconductors [1][2][3][4][5][6]. Low frequency, radiofrequency, microwave and optical reflectivity performed on crystals and thin films give information about the pairing of the superconducting in these materials. In particular the power law dependence of the London penetration depth in the ab plane(T)-(0)~T n , with n~2-2.8 [6] at T<0.5Tc. Such a behaviour is attributed to pair breaking scattering in the so-called s ± superconductivity [6]. Contradicting results on surface resistance have been reported. A coherence peak was observed in thin films with THz and optical measurements [3,4] but no coherence peak was observed in crystals with microwave experiments [1,2]. Here we report radiofrequency surface measurements on optimally doped BaFe 1.9 Ni 0.1 As 2 and overdoped BaFe 1.88 Ni 0.12 As 2 crystals in the frequency range 20MHz -1.5 GHz with a single coil technique.

Careful analysis of the impedance measurements permits us to extract the real part of the conductivity.

Experiment

The crystal were grown using Fe/Ni-As self flux method, details are given in [7]. The samples are platelikes with the plates being perpendicular to the crystallographic c-axis, they exhibit a multicrystal structure with approximately 100 µm size single crystals which are randomly oriented in the ab plane.

Samples were cleaved from a larger crystal. We selected samples with typical dimensions 1000x1000x100 µm 3 , the smallest dimension is along the c axis. The platelet samples are placed inside a copper coil (11 mm length, 2.6 mm diameter, 18 turns, inductance 0.2 µHenry). The coil is situated at the end of a coaxial line inside a terminal adapter. Radio frequency magnetic field is applied parallel to the ab plane. Incident radiofrequency power was fixed to -20dBm. Non resonant measurements of the real (R) and imaginary (L) of the impedance of the coil were performed with an automated impedance analyzer Agilent 4395 in the frequency range 1-100MHz. Measurements of the self LC resonant frequency and series resistance, L=0.2 µH, C=50 fH, quality factor ~ 80, were performed at 1.5 GHz with a Hewlett Packard 8720B network analyzer. Resistance and inductance were measured separately in the absence of a sample and were subtracted from measurements with the sample present [8].

The formulation of the impedance of the coil surrounding the sample was obtained using the equivalent circuit based on a transformer analogy developed in [9]. In this model the primary of the transformer is the measuring coil L 0 . The secondary is defined by an inductance L 2 which is related to the eddy currents induced in the sample, L 2 is a geometrical factor and does not depend on the sample properties. The mutual inductance M between the sample and the coil is defined by the mutual inductance between the primary and secondary, M=k 2 L 0 L 2 , where k is the geometrical coupling factor between the primary and secondary. The inductance of the coil is given by  

2 22 0 2 2 02 0 0 2 2 2 2 ( ) ( ) k L L X L k L L R Z R j L R j X L R j X L                       (1)
Where R and X are the real and imaginary parts respectively of the surface impedance of the sample, R 0 and L 0 are the resistance and inductance of the empty coil and  is the angular frequency.

In the experiments the same coil L 0 =0.2µH was used. Parameters k 2 ~0.1 and L 2 ~0.6 nH were evaluated from the measurements in the normal state of the samples at a temperature T above superconducting transition Tc where the following equality is

0 /2 N N dc XR     (2)
where  0 is the magnetic permeability of vacuum and  dc = 10 6 m the dc sample resistivity measured in [7]. In the superconducting state R and X are deduced from equation (1).

The accuracy of our technique was tested by measurements on superconductor MgB2 powder. Data averaging techniques were used.

Results and Discussion

The temperature dependence of R and X normalized to the value R N measured at 25K are shown for the optimal doped BaFe 1.9 Ni 0.1 As 2 and overdoped BaFe 1.88 Ni 0.12 As 2 crystals in Figs. 1 and2. Below Tc the frequency dependence of X/R N observed is  0.5 which is expected in the superconducting state.

Reactance X is proportional to the London penetration depth and frequency [START_REF] Trunin | Spectroscopy of High Tc superconductors[END_REF], X=   in the superconducting state and R N ~0.5 in the normal state.

The changes in the London penetration depth (T)-(0) deduced from X/R N are shown in Fig. 3.

 follows a ~T2 temperature dependence at T<Tc/2 for optimally doped crystal.

At very low temperatures T<<Tc,  following a power law (T/T C ) n with n~2-2.8 attributed to the effect of strong pair breaking scattering has been reported [6].

After averaging the data, a drop of three orders of magnitude in R/R N is obtained below the superconducting transition Tc. R/R N is proportional to   it results that the expected variation R~   is verified in the superconducting state [START_REF] Trunin | Spectroscopy of High Tc superconductors[END_REF]. Nevertheless the resistance values at T<<Tc in Figs. 1and

2 are larger by more than one order of magnitude than those reported in [1] and [2]. Similar residual surface resistance values are roughly extrapolated for the two different doped crystals in Figs. 1 and3.

Residual loss is a long standing problem in microwave studies of superconductors [START_REF] Bonn | [END_REF]. Residual losses are generally subtracted to obtain the intrinsic resistance. Our resistance measurements are well resolved in the vicinity of the superconducting transition. But at T<<Tc our experimental resolution is not sufficient to extract the intrinsic resistance. The real   and imaginary   part of the conductivity are related to the surface impedance Z by the following relation [START_REF] Trunin | Spectroscopy of High Tc superconductors[END_REF] 0

12 j Z R jX       (3)
The real and imaginary parts of the conductivity normalized to the value      at 25K in the normal state are obtained from the measurements of X/R N and R/R N using the following equations which are deduced from equation (3):

2 1 2 22 / 4 ( / ) ( / ) N N NN R X R R R X R         22 2 2 22 ( / ) ( / ) 2 ( / ) ( / ) NN N NN X R R R R R X R          (4)
  / N and   /(0) , with     T=0), are shown in Figs. 4, 5 and 6.

The conductivity  2 / 2 (0) is related to the London penetration depth [2,[START_REF] Trunin | Spectroscopy of High Tc superconductors[END_REF] and follows a T 2 temperature dependence:

  has been calculated after elimination of the residual loss. The main effect of the subtraction of the residual loss is to force the conductivity to approach zero at zero temperature [START_REF] Bonn | [END_REF]. Regardless of whether or not a residual loss is subtracted, the conductivity   increases with decreasing temperature below Tc to a maximum of about 2-10 times its normal value and this maximum decreases with increasing frequency from 20 to 100MHz.

At the lowest temperatures R/R N is small, it results that the experimental error on   is important.

  depends on the subtracted residual loss and definitive conclusion about the temperature behaviour of   at the lowest temperatures remains not well known.

Our results showing an increasing of   below Tc are in agreement with the measurements on Cobalt doped BaFe 2-x Co x As 2 crystals at 40 GHz which point out a similar behaviour [1]. The authors in [1] argue that the strong rise in the conductivity   below Tc results from a temperature dependent scattering rate decreasing rapidly with temperature.

The phenomenological two fluid model developed in [1,[START_REF] Trunin | Spectroscopy of High Tc superconductors[END_REF] is used here to extract the quasiparticle scattering time . In this two fluid model the temperature dependence of the superconducting density n s is related to the London penetration, n s =    .The normal fluid density n n is given by n n =1-n s , the total carrier density being 1. It results that the temperature dependence of n n in our experiment follows a T 2 dependence, n n ~1-     C ) 2 .

Following the two fluid model, at low frequency   1 is roughly proportional to the product n n Drude like conductivity) which gives in our case

  n n  C ) 2 
Consequentlyscattering rate   is related to the normalized conductivity    N by 2 1 11

N NC T T          (6)
where    is the scattering rate in the normal state at Tc. The scattering rate deduced from    N at 100MHz and 1.5GHz using equation ( 6) is shown in Fig6.

The temperature decrease ofthe scattering rate in the BaFe 2-x Ni x As 2 can be compared with that observed in the Co doped crystals [1]. Below Tc a similar temperature decrease of   is observed in both Ni and Co doped BaFeAs2 superconductors.

Conclusion

Radiofrequency measurements of real and imaginary parts of the conductivity were performed on optimally and overdoped BaFe 2-x Ni x As 2 using a single coil technique. The temperature dependence of the London penetration depth is in agreement with the published results. The increase of the conductivity   below Tc is attributed to the temperature dependence of the quasiparticle scattering rate. A similar behaviour was observed in Co doped BaFe 2-x Co x As 2 .Strong temperature dependence of the quasi particle scattering time is observed in many unconventional superconductors [1,[START_REF] Bonn | [END_REF].

Figure 5

Temperature dependence of the real    N (filled symbols) and imaginary    2 (0) (open symbols) parts of the conductivity normalized to the value  N at 25 K for the over doped BaFe 1.88 Ni 0.12 As 2 crystals. Solid line is calculated with 1-(T/16) 2 .

Figure 6

The temperature dependence of the quasiparticle scattering rate deduced from the conductivity  / N with Eq. 6.

Figure Captions Figure 1 Figure 2

 12 Figure Captions

Figure 3 Figure 4

 34 Figure 3 Temperature dependence of the London penetration depth variation at different frequencies. Figure 4 Temperature dependence of the real    N (filled symbols) and imaginary    2 (0) (open symbols) parts of the conductivity normalized to the value  N at 25 K for the optimally doped BaFe 1.9 Ni 0.1 As 2 crystals. Solid line is calculated with 1-(T/20) 2 .

  Fig1


Aisling Ní Annaidh

Karine Bruyère

Michel Destrade

Michael D Gilchrist

Corrado Maurini

Melanie Otténio

Giuseppe Saccomandi

Automated Estimation of Collagen Fibre Dispersion in the Dermis and its Contribution to the Anisotropic Behaviour of Skin

Keywords: Fibre Orientation, Anisotropic, Skin, Collagen Fibres

Collagen fibres play an important role in the mechanical behaviour of many soft tissues. Modelling of such tissues now often incorporates a collagen fibre distribution. However, the availability of accurate structural data has so far lagged behind the progress of anisotropic constitutive modelling. Here, an automated process is developed to identify the orientation of collagen fibres using inexpensive and relatively simple techniques.

The method uses established histological techniques and an algorithm implemented in the MATLAB image processing toolbox. It takes an average of 15 seconds to evaluate one image, compared to several hours if assessed visually. The technique was applied to histological sections of human skin with different Langer line orientations and a definite correlation between the orientation of Langer lines and the preferred orientation of collagen fibres in the dermis (P<0.001, R 2 = 0.95) was observed. The structural parameters of the Gasser-Ogden-Holzapfel (GOH) model were all successfully evaluated. The mean dispersion factor for the dermis was κ = 0.1404 ± 0.0028. The constitutive parameters µ, k 1 and k 2 were evaluated through physically-based, least squares curve-fitting of experimental test data. The values found for µ, k 1 and k 2 were 0.2014 MPa, 243.6 and 0.1327, respectively. Finally, the above model was implemented in ABAQUS/ Standard and a finite element (FE) computation was performed of uniaxial extension tests on human skin. It is expected that the results of this study will assist those wishing to model skin, and that the algorithm described will be of benefit to those who wish to evaluate the collagen dispersion of other soft tissues.

Introduction

Collagen fibres govern many of the mechanical properties of soft tissues, in particular their anisotropic behaviour. Due to the complex nature of fibre arrangements, such tissues are often represented as either isotropic or transversely isotropic [START_REF] Evans | On the implementation of a wrinkling, hyperelastic membrane model for skin and other materials[END_REF][START_REF] Bischoff | Finite element modeling of human skin using an isotropic, nonlinear elastic constitutive model[END_REF][START_REF] Bischoff | A rheological network model for the contin-uum anisotropic and viscoelastic behaviour of soft tissue[END_REF]. The availability of quantitative structural data on the orientation and concentration of collagen fibres is crucial in order to describe the behaviour of these soft tissues accurately.

A number of studies have used statistical distributions to describe the fibre arrangement in soft tissues. Lanir [START_REF] Lanir | Constitutive equations for fibrous connective tissues[END_REF] was the first to attempt to account for fibre dispersion. Lanir's method expresses the mechanical response in terms of angular integrals. This technique accounts for the contribution of infinitesimal fractions of collagen fibres in a particular orientation. This approach leads to accurate results but it is not a practical option for efficient numerical implementation [START_REF] Cortes | Characterizing the mechanical contribution of fiber angular distribution in connective tissue: Comparison of two modeling approaches[END_REF]. The second approach uses generalised structure tensors (GST), which are assumed to represent the three-dimensional distribution of the fibres. The strain energy function is then calculated by using the average stretch rather than by using a stretch for each individual fibre.

GST is a simpler method than Lanir's, and is easily implemented in FE algorithms [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF].

In this paper we implement the GST method proposed by Gasser et al [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF], described in detail in Section 2.3. In that model, structural parameters such as the fibre dispersion parameter and mean orientation of fibres are required. In order to quantify the orientations of collagen fibres in the human dermis an automated process is developed for detecting collagen orientations in histology slides.

Histology is the microscopic study of cells and tissue. It is an important diagnostic tool in medicine which is utilised here for the purpose of analysing collagen orientation. Traditionally, histology slides are examined by expert observers who assess individual slides visually. This method is both time consuming and subjective. An automated process would enable large volumes of images to be analysed quickly and improve the objectivity of the task. Much research has been conducted recently on imaging techniques of biological soft tissue for the extraction of structural data; however these techniques require either expensive equipment or manual post processing [START_REF] Flamini | Imaging arterial fibres using diffusion tensor imaging -feasability study and preliminary results[END_REF][START_REF] Verhaegen | A reliable, non-invasive measurement tool for anisotropy in normal skin and scar tissue[END_REF][START_REF] Wu | Automated quantification and reconstruction of collagen matrix from 3D confocal datasets[END_REF][START_REF] Yasui | Characterization of collagen orientation in human dermis by two-dimensional second-harmonic-generation polarimetry[END_REF]. There has been little research to date carried out on the automated analysis of histology slides. Four notable exceptions are the work of Van Zuijlen et al [START_REF] Van Zuijlen | Morphometry of dermal collagen orientation by Fourier analysis is superior to multi-observer assessment[END_REF], Noorlander et al [START_REF] Noorlander | A quantitative method to determine the orientation of collagen fibers in the dermis[END_REF], Elbischger et al [START_REF] Elbischger | Automatic analysis of collagen fiber orientation in the outermost layer of human arteries[END_REF] and Jor et al [START_REF] Jor | Modelling collagen fibre orientation in porcine skin based upon confocal laser scanning microscopy[END_REF].

Elbischger et al [START_REF] Elbischger | Automatic analysis of collagen fiber orientation in the outermost layer of human arteries[END_REF] successfully automated the analysis of collagen fibre orientation in the human adventitia, (the outermost membrane of the artery). Their technique was based on the automated analysis of Transmitted Light Microscopy images which were stained with Elastica van Gieson. The algorithm uses a ridge and valley analysis to detect the orientation of collagen fibres and segments regions of homogeneous fibre orientations. This is a somewhat complicated technique which requires substantial coding to implement. There are also assumptions made which cannot be applied to the dermis, i.e. that the image contains at least 50% collagen and that the fibres have a common orientation. In contrast, Noorlander et al [START_REF] Noorlander | A quantitative method to determine the orientation of collagen fibers in the dermis[END_REF] used a relatively simple technique upon which our own algorithm is based. In their study, Picrosirius red staining was used together with epipolarised light to image the fibres. Individual collagen fibres were detected and ellipses fit to the longest ten fibres. While their technique is quantitative it is not an automated technique and requires the user to identify fibre orientations visually. The results refer to the 'collagen alignment index' which is the mean length of the best fit ellipse and has no connection with the angular orientation of fibres, which itself is not measured. Van Zuijlen et al [START_REF] Van Zuijlen | Morphometry of dermal collagen orientation by Fourier analysis is superior to multi-observer assessment[END_REF] used Fourier analysis to measure the level of anisotropy of collagen in the skin.

Those authors concluded that analysis of 'orientation index' by Fourier analysis is superior to conventional techniques. However, the 'orientation index' is a measure of the anisotropy of the matrix but, just like the 'collagen alignment index' of Noorlander et al [START_REF] Noorlander | A quantitative method to determine the orientation of collagen fibers in the dermis[END_REF], it fails to provide information on the mean orientation of the fibres. In a recent publication Jor et al [START_REF] Jor | Modelling collagen fibre orientation in porcine skin based upon confocal laser scanning microscopy[END_REF] modelled the orientation of collagen fibres in porcine skin. Their technique involved the staining of porcine skin sections with Picrosirius red and image acquisition by confocal laser scanning microscopy. Their study provided quantitative results on the orientation of collagen fibres in porcine skin, however the plane of interest was normal to the epidermis and it is generally believed that the preferred orientation of collagen fibres is parallel to the epidermis [START_REF] Holzapfel | Handbook of Materials Behavior Models: Biomechanics of Soft Tissue[END_REF].

The chief advantage of the technique proposed here over other automated techniques is that it is an inexpensive and relatively simple technique. It can be easily implemented in MATLAB using the Image Processing Toolbox so that the algorithm can be easily amended as required for the user's specific application. It is a fully automated process (aside from the image acquisition phase) and is capable of analysing a greater number of images, and faster, than a manual analysis. It also eliminates the subjectivity which is present using traditional methods. The technique provides quantitative data on both the collagen orientation and the level of anisotropy in the dermis, but it could be easily amended to identify any tissues/objects within other soft tissues.

Due to the anisotropic properties of collagen based soft tissues [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF], modelling now often incorporates a collagen fibre distribution. The accuracy of these structural models rely heavily on knowledge of collagen fibre orientations. While this data has recently been published for porcine skin by Jor et al [START_REF] Jor | Modelling collagen fibre orientation in porcine skin based upon confocal laser scanning microscopy[END_REF], to the best of the authors' knowledge, quantitative data on the orientation of collagen fibres in the human dermis has not been previously published. With the quantitative data obtained here, structural parameters such as γ and κ from the Gasser-Ogden-Holzapfel model [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF] can be evaluated independently. This structural data, coupled with the mechanical testing of the same skin samples [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF] provides us with sufficient data to model human skin using the GOH model.

Materials & Methods

Tensile tests and histology

In vitro tensile tests of human skin were performed using a Universal Tensile Test machine at a strain rate of 0.012s -1 . The tensile load was measured with a 1kN piezoelectric load cell and the strain was measured via a displacement actuator. Skin biopsies were excised adjacent to the tensile test samples and stained with Van Gieson to dye collagen fibres pink/red. Further details of the tensile test experimental procedure and histology protocol are outlined in Ní Annaidh et al [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF].

Automated detection of fibre orientations

Images were taken of slides parallel to the epidermis (see Section 2.4) using an Aperop ScanScope XT scanner and ScanScope software. This slide scanner has the ability to scan up to 120 slides automatically. The images to be analysed were selected from the reticular layer of the dermis which forms the main structural body of the skin. The optical magnification used was 5x, which is quite low; however, this was the optimal magnification to obtain a large enough field of view to be representative while also maintaining enough detail to recognise the boundaries of fibre bundles. Multiple images were taken of each sample in order to capture the entire area.

The field of view for each image taken at this magnification was 2.24 mm 2 . The orientation of collagen fibres were then calculated in a fully-automated customised MATLAB routine using the Image Processing Toolbox. The algorithm is described below and is illustrated in Fig. 1.

• A global threshold level was computed using the graythresh function. This function chooses a global image threshold using Otsu's method which aims to minimise the intraclass variance of black and white pixels [START_REF] Otsu | A threshold selection method from grey-level histograms[END_REF]. This automatically distinguishes collagen fibres from other areas which are not of interest in the analysis, such as cells or histological ground substance. The image was then binarised based on this threshold level producing an image containing black pixels for collagen and white pixels for all other areas;

• Morphological operations were performed on the binary image using the bwmorph function. An erosion step was performed to detach cross-linking fibres from each other. This step removed pixels from the boundaries using a structural element of size 3 pixels x 3 pixels. One iteration only of this step was performed which was the optimal number of iterations to detach cross-linking fibres while also leaving smaller fibres intact. The second morphological operation was performed using the imfill function. This function fills in the 'holes' within the binary image, where a hole is defined as a set of isolated pixels which cannot be reached by filling in the background pixels from the edge. This step resulted in a 'cleaner' image to analyse.

• Individual fibre bundles are identified using bwlabel. This function identifies connected components (8-connected) in a 2D binary image and labels each component individually.

It uses the general procedure outlined by Haralick and Shapiro [START_REF] Haralick | Computer and robot vision[END_REF];

• The regionprops function outputs a set of properties for each labelled component. The function also fits an ellipse to each component by matching the second order moments of that component to an equivalent ellipse following the procedure by Haralick and Shapiro [START_REF] Haralick | Computer and robot vision[END_REF];

• Only ellipses which were elongated and of a certain area were selected 1 ;

• The orientation of the major axis of each ellipse was calculated and taken as the approximate orientation of each component. The advantage of fitting an ellipse about each component is that the orientation is measured in a systematic and repeatable manner which can account for the non-uniformities of the shape of each component;

• The orientations of each component was then plotted on a histogram (see Fig. 2). Two distinct peaks were evident from this figure. It is assumed that these two peaks correspond to the preferred orientation of two crossing families of fibres as shown in Fig. 3. A von Mises probability density function was fit to the data and the mean orientation and dispersion factor were calculated as described in Section 2.4.

Validation

As explained by Van Zuijlen et al [START_REF] Van Zuijlen | Morphometry of dermal collagen orientation by Fourier analysis is superior to multi-observer assessment[END_REF], there are difficulties involved with validating automated techniques because the true result of each histology slide is unknown. For the validation of their algorithm, computer generated images representing collagen fibres with known orientations were used. However, in order to create an accurate computer generated image, a model which represents the structure of the collagen matrix is required, but no such model exists [START_REF] Elbischger | Automatic analysis of collagen fiber orientation in the outermost layer of human arteries[END_REF]. To validate our technique, a selection of slides in the plane perpendicular to the epidermis (see Fig. 4) were manually segmented and their mean orientation compared to those calculated through the automated process. Although it is the plane parallel to the epidermis that is used for the collection of structural data, the perpendicular plane was chosen for validation purposes because in general it displays a more orientated pattern and is easier to segment manually. 1 Out-of-plane fibres manifest themselves as circular areas. This condition was introduced to exclude out-ofplane fibres from the analysis. A sensitivity analysis was carried out to determine the sensitivity of the algorithm to varying the area and eccentricity criteria. It was found that a ±10% change in both the area and eccentricity criteria led to a ±1% change in the mean orientation. This indicates that the selected criteria do not significantly affect the result of mean orientation. For this study it was specified that in-plane fibres must have an eccentricity larger than 0.7 and an area greater than 1000 pixels, which for our images, captured at 5x, corresponds to 5µm 2 .

Manual segmentation was performed by marking bundles of elongated fibres with a line as shown in Fig. 5. The orientation of each line was then measured manually and the mean calculated.

The dermis possesses a complex interwoven collagen structure, which is difficult to assess.

It is quite possible that the collagen pattern may change over a small volume. For this reason six images were taken at each depth, whereby each image spanned 3.75mm 2 of the 50mm 2 sample. The analysis was performed at six levels of increasing depth and the mean orientation was taken as the average over these six levels.

Gasser-Ogden-Holzapfel Model

The Gasser-Ogden-Holzapfel (GOH) model applies to incompressible solids with two preferred directions aligned along the unit vectors a 1 and a 2 (say) in the reference configuration (see Fig. 3). Its strain energy density Ψ is of the form

Ψ = Ψ(C, H 1 , H 2 ) ( 1
)
where C is the right Cauchy-Green strain tensor, and the structure tensors H 1 , H 2 depend on a 1 and a 2 and on the dispersion factors κ 1 , κ 2 (to be detailed later), respectively, as follows

H i = κ i I + (1 -3κ i )a 1 ⊗a 2 , (i = 1, 2). (2)
Specifically, the GOH model assumes that Ψ depends on I 1 = tr(C), tr(H 1 C) and tr(H 2 C) only, as follows

Ψ = µ 2 (I 1 -3) + µ ∑ i=1,2 k i1 2k i2 e k i2 [tr(H i C)-1] 2 -1 , (3)
where µ, k i1 , k i2 are positive material constants, and from Eq. 2, tr(

H i C) = κ i I i + (1 -3κ i )I 4i , (4)
with I 4i ≡ a i •Ca i , two anisotropic invariants. Note that the constitutive parameter µ has the dimensions of stress: it would be the shear modulus of the solid if there were no fibres (k i1 = 0); whilst the parameters k i1 and k i2 are dimensionless stiffness parameters: the k i1 are related to the relative stiffness of the fibres in the small strain regime, and the k i2 are related to the large strain stiffening behaviour of the fibers. Now we focus on homogeneous uniaxial tensile tests. These can be achieved for anisotropic tissues when two families of fibres are mechanically equivalent k 11 = k 21 ≡ k 1 (say) and k 12 = k 22 ≡ k 2 (say), with the same dispersion factor κ 1 = κ 2 ≡ κ (say), and when the tension occurs along the bisector of a 1 and a 2 (see Fig. 3). Let us call γ the angle between a 1 and the tensile direction, so that now

a 1 = cos γ i + sin γ j, a 2 = cos γ i -sin γ j. (5)
Here, i is the unit vector in the direction of tension, and j is the unit vector in the lateral direction, in the plane of the sample. The stretch ratios along those unit vectors are λ 1 and λ 2 , respectively. Then I 41 = I 42 = λ 2 1 cos 2 γ + λ 2 2 sin 2 γ ≡ I 4 (say) and Ψ reduces to

Ψ = µ 2 (I 1 -3) + µ k 1 k 2 e k 2 [κI 1 +(1-3κ)I 4 -1] 2 -1 , (6)
giving the following expression for σ , the Cauchy stress tensor

σ = -pI + 2 ∂ Ψ ∂ I 1 FF T + ∂ Ψ ∂ I 4 [Fa 1 ⊗Fa 1 + Fa 2 ⊗Fa 2 ] , (7)
where p is a Lagrange multiplier introduced by the internal constraint of incompressibility and

F is the deformation gradient. Note that Fa 1 ⊗Fa 1 +Fa 2 ⊗Fa 2 = 2(λ 1 cos γ) 2 i ⊗ i+2(λ 2 sin γ) 2 j ⊗ j,
showing that σ is diagonal in the {i, j, k} basis. Its components are

σ 11 = -p + 2(Ψ 1 + Ψ 4 cos 2 γ)λ 2 1 = 0, σ 22 = -p + 2(Ψ 1 + Ψ 4 sin 2 γ)λ 2 2 = 0, σ 33 = -p + 2Ψ 1 λ -2 1 λ -2 2 = 0, (8)
where

2Ψ 1 = µ(1 + 4k 1 καe k 2 α 2 ), 2Ψ 4 = 4µk 1 (1 -3κ)αe k 2 α 2 , α = κ(λ 2 1 + λ 2 2 + λ -2 1 λ -2 2 ) + (1 -3κ)(λ 2 1 cos 2 γ + λ 2 2 sin 2 γ) -1. (9)
Now, eliminate p from the stress components to get the two equations

σ 11 = µ(λ 2 1 -λ -2 1 λ -2 2 ) + 4µk 1 αe k 2 α 2 κ(λ 2 1 -λ -2 1 λ -2 2 ) + (1 -3κ)λ 2 1 cos 2 γ , (10)
0 = λ 2 2 -λ -2 1 λ -2 2 + 4k 1 αe k 2 α 2 κ(λ 2 2 -λ -2 1 λ -2 2 ) + (1 -3κ)λ 2 2 sin 2 γ . (11)
Equation [START_REF] Gamage | Identification of mechanical properties of heterogeneous soft bodies using gravity loading[END_REF] gives the relationship between the tensile stretch and the lateral stretch and allows, implicitly, λ 2 to be expressed in terms of λ 1 . Substituting then into [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF] gives the σ 11λ 1 stress-stretch relationship. In the isotropic limit, κ = 1/3 (see Section 2.4), and ( 11) yields the well-known relationship

λ 2 = λ -1/2 1
for uniaxial tension in incompressible solids.

The two equations ( 10)- [START_REF] Gamage | Identification of mechanical properties of heterogeneous soft bodies using gravity loading[END_REF] form the basis of a numerical determination of the constitutive parameters µ, k 1 and k 2 , assuming that the structural parameters κ and γ are known. It should be noted here that the inclusion of µ in the anisotropic term of equation ( 6) is not standard for the GOH model and has been added here for ease of calculation. We now quantify further those latter parameters, κ and γ .

Fibre Dispersion

The GOH model assumes that the mean orientation of collagen fibres has no out-of-plane component. Our histological examination of the skin indicates that the majority of collagen fibres in the dermis run parallel to the epidermis. Slides parallel to the epidermis had, on average, three times less cross-sectioned fibres than slides perpendicular to the epidermis. This is in agreement with Holzapfel et al [START_REF] Holzapfel | Handbook of Materials Behavior Models: Biomechanics of Soft Tissue[END_REF] who state that the preferred orientation of the 3D collagen fiber network lies parallel to the surface, but to prevent out-of-plane shearing, some fiber orientations have components which are out-of-plane. Despite the assumption that the fibres have no out-of-plane component in the GOH model, the three-dimensional nature of the adopted distribution implies that although the preferred orientation of the fibers are in the plane parallel to the epidermis, some fibers orientations have an out-of-plane component [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF].

Here we assumed that each of the two families of collagen fibres is distributed according to a π-periodic Von Mises distribution, which is commonly assumed for directional data. The standard π-periodic Von Mises Distribution is normalized and the resulting density function, ρ(Θ), reads as follows,

ρ(Θ) = 4 b 2π exp[b(cos(2Θ) + 1] erfi( √ 2b) , ( 12
)
where b is the concentration parameter associated with the Von Mises distribution and Θ is the mean orientation of fibres (for a graph of the variation of the dispersion parameter κ with the concentration parameter b, see Gasser et al [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF]).

The parameters b and θ were evaluated using the mle function in MATLAB. Analogous to least squares curve-fitting, the maximum likelihood estimates (MLE), is the preferred technique for parameter estimation in statistics. Then κ is calculated by numerical integration of the integral given by Gasser et al [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF],

κ = 1 4 π 0 ρ(Θ) sin 3 ΘdΘ. ( 13
)
The structural parameter κ describes the material's degree of anisotropy. It must be in the range 0 κ 1/3: the lower limit, κ = 0, relates to the ideal alignment of collagen fibres and the upper limit, κ = 1/3, relates to the isotropic distribution of collagen fibres. Fig. 6 is a 3D graphical representation of the orientation of collagen fibres for different values of κ.

The fibres were assumed to form an interweaving lattice structure as first postulated by Ridge and Wright [START_REF] Ridge | Mechanical properties of skin: A bioengineering study of skin structure[END_REF] and shown in Fig. 3. These authors suggested that the mean angle of the two families of fibres indicates the direction of the Langer lines. More recent in vitro [START_REF] Jor | Modelling collagen fibre orientation in porcine skin based upon confocal laser scanning microscopy[END_REF] and in vivo [START_REF] Ruvolo | Skin Viscoelasticity Displays Site-and Age-Dependent Angular Anisotropy[END_REF] studies have also supported this hypothesis. The lattice structure proposed by Ridge and Wright [START_REF] Ridge | Mechanical properties of skin: A bioengineering study of skin structure[END_REF] is an idealised one, and the adoption here of the dispersion factor creates a more realistic scenario. Finally, recall that the two families (directions) of fibres are assumed to have a common dispersion factor.

Finite element representation

An FE computation was used to simulate uniaxial tensile tests of human skin which were described in a previous publication [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF]. The simulation was carried out for three samples: parallel, perpendicular, and at 45 • to the Langer lines. The test samples were of the dimensions shown in Fig. 7. The length of the unclamped specimen was 68 mm and the thickness was 2.25 mm. 1512 reduced integration hybrid hexahedral (C3D8RH) elements were used for the mesh. The numerical analyses were performed using the static analysis procedure in ABAQUS/Standard. Displacement was applied through a smooth amplitude boundary condition, and the top and bottom of the sample were encastred to represent the clamping of samples.

The material model used was the anisotropic GOH model, which is an internal material model in ABAQUS.

3 Results

Structural parameters

As expected, it was found from the histology that the collagen fibres were locally orientated.

This meant that each sample had a different mean orientation and fibre dispersion. Table 1 tabulates the results of 12 different human skin samples, with different orientations, which were procured from the backs of two different subjects2 . See Fig. 8 for details of specimen location and orientation. The preferred orientation, Θ, refers to the bisector of the two families of fibres. The 95% Confidence Interval of the mean of these images ranged from 1.16 • -2.77 • , thereby indicating that the preferred orientation does not change significantly over this small area. The results of our validation revealed that the automated process differed from the manual segmentation by an average of 5 • ±4 • , giving us further confidence in the validity of our technique.

A Pearson correlation test was carried out to test for a correlation between the measured preferred orientation obtained through histology and the perceived orientation of Langer lines (the natural lines of tension in the skin). The orientation of Langer lines was assessed using generic maps, described further in Ní Annaidh et al [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF].The correlation was deemed to be significant (P<0.001) with an R 2 value of 0.95. This shows that the Langer lines have an anatomical basis, a point which had previously been suggested but until now had not been quantitatively assessed.

Constitutive parameters

It was assumed that the orientation of collagen fibres is symmetric about the axis of applied stress. In reality, this is not always the case; however some assumptions must be made in order to ensure that the constitutive relations remain practical for numerical implementation.

In particular, as explained in Section 2.3, this assumption leads to a homogeneous deformation of the sample, and in turn, to an explicit stress-strain solution. In this section, three illustrative examples from Table 1 (highlighted) have been chosen for further investigation. Two of these examples have been chosen because they are the samples that are closest to being symmetrical about the axis of applied stress. The third sample has been chosen to illustrate how its behaviour can be modelled using FE analysis and is described further in Section 3.3

The constitutive parameters for the GOH model are obtained by using Equations ( 10) and [START_REF] Gamage | Identification of mechanical properties of heterogeneous soft bodies using gravity loading[END_REF], obtained in Section 2.3. When linearized in the neighbourhood of small strains, λ i 1 + ε i , say, we find that they read as follows,

σ 11 = 4µ[1 + 2k 1 (1 -3κ) 2 cos 4 γ]ε 1 + 2µ[1 + 4k 1 (1 -3κ) 2 sin 2 γ cos 2 γ]ε 2 , ( 14
) 0 = [1 + 4k 1 (1 -3κ) 2 cos 2 γ sin 2 γ]ε 1 + 2[1 + 2k 1 (1 -3κ) 2 sin 4 γ]ε 2 . ( 15
)
These expressions reveal that the constitutive parameters µ and k 1 are related to the early stages of the tensile tests, whilst k 2 is a stiffening parameter, related to the latter (nonlinear) stages of the tensile tests. By solving (15) for ε 2 , and substitution into [START_REF] Holzapfel | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF], we find the linear stress-strain relation σ 11 = E 1 ε 1 , where E 1 is the infinitesimal Young modulus in the 1-direction, found here as

E 1 = 3 + 8k 1 (1 -3κ) 2 (1 -3 cos 2 γ sin 2 γ) 1 + 2k 1 (1 -3κ) 2 sin 4 γ µ, (16)
(which is consistent with the formula E = 3µ in linear isotropic (κ = 1/3) incompressible elasticity). Hence, by plotting the values of σ 11 for the early part of the tests (first 1000 data say, corresponding to a tensile stretch of less than 2%), we can determine E 1 by linear regression analysis, see Fig. 9(a). Here we have plotted the 'parallel' sample highlighted in Table 1 data for which was collected from tensile tests of human skin samples [START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF].

Once E 1 is determined, µ can be expressed in terms of E 1 and k 1 using Eq. ( 16). Then the remaining material parameters k 1 and k 2 are found through the nonlinear least squares fitting with experimental test data of Equation ( 10), subject to the definition of λ 2 in terms of λ 1 given by Equation [START_REF] Gamage | Identification of mechanical properties of heterogeneous soft bodies using gravity loading[END_REF]. The data fitting was performed using the lsqnonlin MATLAB routine in the Optimisation Toolbox where the objective function, Err(k), was given as

Err(k) = n ∑ i=1 (y exp i -y model(k) i ) 2 ( 17
)
Where n is the number of experimental data points, y exp i is the experimental value and

y model(k) i
is the value predicted by the model using the current material parameters, k.

Non-linear optimisation procedures are often sensitive to the initial starting point provided by the user [START_REF] Ogden | Fitting hyperelastic models to experimental data[END_REF]. In our case, the initial estimate for k 1 was found by calculating the slope of the non-linear part of the stress-stretch curve. Since we have shown that k 1 is related to the stiffening stage of the tensile test, our initial estimate therefore has a physical meaning.

Furthermore, the initial estimates of both k 1 and k 2 were varied over a large range and lead to the same set of optimal parameters each time, illustrating that the results of the optimisation procedure are not sensitive to this initial estimate.

The results of our optimisation procedure gave us a value of 243.6 for k 1 and 0.1327 for k 2 , with an R 2 of 99.5%. Fig. 9(b) shows the GOH model fit to the experimental data. It can be seen that these material parameters provide an excellent fitting to the 'parallel' sample, at least from a descriptive point of view. Turning now to the predictive capabilities of the GOH model, we examine a sample 'perpendicular' to the Langer lines. In Fig. 10 we use the material parameters k 1 and k 2 obtained through least squares fitting, µ calculated by linear regression and Eq. ( 16), coupled with the unique structural data for the 'perpendicular' sample in Table 2.

We compare the model prediction for a tensile test occurring perpendicular to the Langer lines to the experimental data. The fit remains good with R 2 = 97.96%. This shows that the model is capable of predicting the behaviour of skin once the structural parameters have been evaluated.

Finite element simulations

The conventions used by ABAQUS are related to ours through

C 10 = µ/2, k 1 = k 1 µ and k 2 = k 2 ,
so that the ABAQUS parameters used were C 10 = 0.1007 MPa, k 1 = 24.53 MPa and k 2 = 0.1327.

The FE simulation results of the uniaxial tensile tests for both the parallel and perpendicular samples were identical to the analytical solution. The results were independent of both mesh density and element type. Fig. 11 shows the Cauchy stress distribution across the parallel and perpendicular samples at the end of the test. The large difference in magnitudes between the two is due to the variation in the mean orientation of fibres. Note the uniform distribution of stress in the middle section of the test samples, thanks to the dog-bone shape of the specimen, and the symmetry of fibres about the axis of applied stress.

As discussed in Section 3.2, for ease of determining the material parameters, it was assumed that the orientation of collagen fibres is symmetric about the axis of applied stress. This is for an idealised scenario only, where one knows the exact orientation of collagen fibres prior to testing, and can therefore apply the stress in this orientation. However, with the model parameters that have now been determined, a FE simulation can provide results for samples where the collagen fibres are not symmetric about the axis of applied stress. Examining the nonsymmetric example in Fig. 12(a), we can see a non-uniform distribution of stress throughout the test specimen. A local magnification of stress occurs near the neck regions of the test specimen. This is due to the non-symmetry of collagen fibres about the axis of applied stress and makes this problem a much more complicated one to solve analytically. The presence of significant levels of shear in Fig. 12(b) (which is absent from both the parallel and perpendicular samples)

indicates further the effect of this non-symmetry on the sample response.

Because the stress distribution in this sample is non-uniform, to examine the predictive capabilities of the GOH model here, we must plot the experimental force-displacement data against the values predicted by ABAQUS for a node at the top of the test sample, see Fig. 13.

Again, we have found that the model predicts the behaviour well with an R 2 of 94.4%, showing that the GOH model is capable of predicting the anisotropic response of human skin.

Discussion

For this paper, histology slides in three different planes were examined; however, after capturing the images from all three planes it was observed that three times as many cross-sectioned fibres were present in the plane normal to the epidermis, therefore an assumption was made to ignore the fibres normal to the epidermis. Hence, the further analysis of the samples was restricted to the plane parallel to the epidermis alone, making this a 2D analysis. Physically however, there are a percentage of fibres that run normal to the epidermis and this information has not been captured here. Furthermore, unloaded collagen fibres have a crimped nature: Their relative orientation may vary with respect to the epidermal plane. Here we have tried to overcome this limitation by taking an average measure over six levels of depth spanning 30µm,with the view that the general orientation of the collagen fibres are still captured. Ideally, a full 3D analysis of the dermal structure could consider the effects of collagen crimping. A 3D analysis was beyond the scope of this paper, but the current technique could be extended by creating a montage of overlapping images through the thickness of the dermis, as described in Jor et al [START_REF] Jor | Modelling collagen fibre orientation in porcine skin based upon confocal laser scanning microscopy[END_REF], therefore turning the 2D analysis into a 3D analysis.

A further limitation of this technique is that we have assumed that the structure of the skin biopsy removed is representative of the entire tensile test sample. In reality, the structure, and therefore the properties of skin may vary considerably over a small area. The variation of the mean fibre orientation over the volume of the biopsy was quantified, however, the size of the biopsy is very small relative to the tensile test sample and we cannot infer that the variation would be negligible. At a minimum, future studies should excise a number of biopsies along the length of the test specimen to investigate the variation of the structure.

While this technique has been described as 'automated', there are still a number of 'manual' steps that must first be performed. The first manual task is the histological staining, combined with the mounting of skin biopsies. The collagen detection process demands a high quality of histological staining for the method to be successful and therefore, care must be taken to follow standard procedures carefully. The second 'manual' task is the image acquisition phase.

Modern 'slide scanners' automatically scan multiple slides at once meaning that tedious image acquisition techniques using a manual microscope are no longer necessary, however images must still be captured from the 'digital slide'.

While this technique has provided quantitative structural data of human skin, it can, of course, only be applied in-vitro. Considering the effect that the mean orientation of collagen fibres has on the mechanical response of skin, the development of in-vivo methods for establishing the orientation of fibres is of the utmost importance. Advanced imaging techniques such as ultrasonic surface wave propagation may eventually provide real-time, in-vivo structural data.

It should be noted that the nonlinear curve fitting technique rests only on measurements of λ 1 and σ 11 , and that λ 2 , along with k 1 and k 2 are obtained during the simultaneous optimisation of Eq. ( 10) and Eq. ( 11). Ideally, a more complete analysis would include, compare, and contrast experimental data for λ 2 and/or λ 3 . An extensive experimental data set would include planar biaxial tests with in-plane shear and separate through thickness shear tests [START_REF] Flynn | Modeling the Mechanical Response of In Vivo Human Skin Under a Rich Set of Deformations[END_REF][15], however in the absence of these advanced testing protocols tensile tests coupled with a histological study of the collagen fibre alignment can be used for reasonable determination of material parameters [START_REF] Holzapfel | Determination of material models for arterial walls from uniaxial extension tests and histological structure[END_REF]. Of course non-uniqueness of optimal material parameters is an intrinsic problem in non-linear fitting. It is possible that the optimisation procedure finds a local minimum and assumes that this is the global minimum [START_REF] Ogden | Fitting hyperelastic models to experimental data[END_REF]. To ensure that the optimisation procedure is providing a unique set of material parameters a number of checks are available: The properties of the Hessian matrix can be investigated at the optimum [START_REF] Gamage | Identification of mechanical properties of heterogeneous soft bodies using gravity loading[END_REF], or alternatively, one can plot the objective function as a function of the varying material parameters. In this case the objective function was investigated and the plot (not reproduced) shows that our procedure calculates a global minimum and not merely a local minimum.

In this study we have developed a simple automated process which can detect the orientation of collagen fibres. This technique can be easily implemented in MATLAB and can be adapted to detect other biological features, such as certain cells, leading to applications in diagnostics.

We have applied this technique to skin biopsies and provided new quantitative data on the ori-entation of collagen fibres in the human dermis. So far, the availability of accurate structural data has lagged behind the progress of anisotropic constitutive modelling. Here we have provided the structural data required to accurately make use of advances in constitutive modelling, and help fill the void of experimental data. The model parameters of the GOH model have been evaluated for skin using experimental data from the same skin samples. These sets of parameters will provide invaluable data for those wishing to model the anisotropic behaviour of skin.

Finally, an FE simulation of a uniaxial tensile test on three separate human skin samples was performed which predicted the the response of these three samples well. We have illustrated that the Gasser-Ogden-Holzapfel model can successfully model the anisotropic behaviour of human skin and that it can be implemented in ABAQUS with ease.

  (a) Original histology slide. Scale bar is 1mm. (b) Binarised image after automated thresholding. (c) Binarised image after erosion step. (d) All identified collagen bundles outlined in green. (e) Remaining bundles which meet area and eccentricity criteria.

  (f) Best fit ellipse about each fibre that meets the specified criteria.

Figure 1 :Figure 2 :

 12 Figure 1: Images output from automated algorithm.

Figure 3 :

 3 Figure 3: Lattice structure of crossing collagen fibres with fibre dispersion taken into account.

Figure 4 :

 4 Figure 4: Biopsies of skin samples for purpose of histological staining. Note that the biopsies have been sliced in three orthogonal planes.

Figure 5 :

 5 Figure 5: Manual segmentation of collagen. Elongated fibres were marked by black lines, and their orientation was later measured manually.

  (a) κ = 0.0085 (b) κ = 0.25 (c) κ = 0.33

Figure 6 :

 6 Figure 6: Three dimensional representation of the orientation of collagen fibres.

Figure 8 :

 8 Figure 8: Location of tensile test samples shown in Table1(figure amended from[START_REF] Langer | On the anatomy and physiology of the skin. The Imperial Academy of Science[END_REF]).

Figure 9 :

 9 Figure 9: Nonlinear curve fitting to obtain the constitutive parameters: µ and k 1 are related to the early (infinitesimal) stress-strain part of the graph, see (a); k 2 , to the rest of the curve.

Figure 10 :

 10 Figure 10: Comparison of sample parallel to Langer lines and perpendicular to Langer lines.

Figure 11 :

 11 Figure 11: Cauchy stress in Pa of sample strained by 50% (a) Parallel to the Langer lines (b) Perpendicular to the Langer lines.

Figure 12 :

 12 Figure 12: Cauchy stress in Pa of a non-symmetric sample strained by 30% (a) σ 11 (b) σ 12

Figure 13 :

 13 Figure 13: Comparison between predicted model response and experimental forcedisplacement data for a non-symmetric sample.

Table 2 :

 2 Values obtained through curve-fitting for the parameters µ, k 1 and k 2 . Also displayed is R 2 , a measure of goodness of fit, and κ and γ obtained directly through histology.

	Sample	µ	k 1	k 2	γ	Orientation Θ	κ	R 2
		MPa				•		%
	Parallel	0.2014 243.6 0.1327	41	88	0.1535 99.54
	Perpendicular	0.2014 243.6 0.1327	41	0	0.1456 97.96
	Non-symmetric 0.2014 243.6 0.1327	41	118	0.1602 94.40

It should be noted here that these results have been obtained using the mle method described in Section 2.4. An alternative, but less reliable method exists whereby data is clustered into two intersecting fibre families using an agglomerative clustering algorithm such as was performed in Ní Annaidh et al[START_REF] Ní Annaidh | Otténio Characterising the anisotropic mechanical properties of excised human skin[END_REF].

Acknowledgements

The authors acknowledge gratefully the advice and assistance of Mr. Ciaran Driver, Dr. Michael Curtis and Prof. Marie Cassidy, of the Office of the State Pathologist (Ireland), in the area of histology. This research was supported by a Marie Curie Intra European Fellowship within the 7 th European Community Framework Programme, awarded to MD; by the Irish Research Council for Science, Engineering and Technology; by the Office of the State Pathologist (Irish Department of Justice and Equality); and by the Ile-de-France region. G.S. is supported by the PRIN 2009 project "Matematica e meccanica dei sistemi biologici e dei tessuti molli".

Tables & Figures