Zheng Cheng
email: zheng.cheng@inria.fr

Massimo Tisi
email: massimo.tisi@inria.fr

A Deductive Approach for Fault Localization in ATL Model Transformations

Keywords: rule RS2RS {, } rule IS2RS {, } rule T2TC {

In model-driven engineering, correct model transformation is essential for reliably producing the artifacts that drive software development. While the correctness of a model transformation can be specified and checked via contracts, debugging unverified contracts imposes a heavy cognitive load on transformation developers. To improve this situation, we present an automatic fault localization approach, based on natural deduction, for the ATL model transformation language. We start by designing sound natural deduction rules for the ATL language. Then, we propose an automated proof strategy that applies the designed deduction rules on the postconditions of the model transformation to generate sub-goals: successfully proving the sub-goals implies the satisfaction of the postconditions. When a subgoal is not verified, we present the user with sliced ATL model transformation and predicates deduced from the postcondition as debugging clues. We provide an automated tool that implements this process. We evaluate its practical applicability using mutation analysis, and identify its limitations.

context HSM!Transition inv Pre1

Introduction

Model-driven engineering (MDE), i.e. software engineering centered on software models and their transformation, is widely recognized as an effective way to manage the complexity of software development. One of the most widely used languages for model transformation (MT) is the AtlanMod Transformation Language (ATL) [START_REF] Jouault | ATL: A model transformation tool[END_REF]. Like several other MT languages, ATL has a relational nature, i.e. its core aspect is a set of so-called matched rules, that describe the mappings between the elements in the source and target model.

With the increasing complexity of ATL MTs (e.g., in automotive industry [START_REF] Selim | Model transformations for migrating legacy models: An industrial case study[END_REF], medical data processing [START_REF] Wagelaar | Using ATL/EMFTVM for import/export of medical data[END_REF], aviation [START_REF] Berry | Synchronous design and verification of critical embedded systems using SCADE and Esterel[END_REF]), it is urgent to develop techniques and tools that prevent incorrect MTs from generating faulty models. The effects of such faulty models could be unpredictably propagated into subsequent MDE steps, e.g. code generation. Therefore, the correctness of ATL is our major concern in this research. Typically correctness is specified by MT developers using contracts [9-13, 19, 21, 23]. Contracts are pre/postconditions that express under which condition the MT is considered to be correct. In the context of MDE, contracts are usually expressed in OCL [START_REF]Object Management Group: The Object Constraint Language Specification[END_REF] for its declarative and logical nature.

In [START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF], Cheng et al. developed the VeriATL verification system to deductively verify the correctness of ATL transformations w.r.t. given contracts. VeriATL automatically generates the axiomatic semantics of a given ATL transformation in the Boogie intermediate verification language [START_REF] Barnett | Boogie: A modular reusable verifier for object-oriented programs[END_REF], combined with a formal encoding of EMF metamodels [START_REF] Steinberg | EMF: Eclipse modeling framework[END_REF] and OCL contracts. The Z3 automatic theorem prover [START_REF] De Moura | Z3: An efficient SMT solver[END_REF] is used by Boogie to verify the correctness of the ATL transformation.

However, when a contract on the MT is not verified, current verification systems like VeriATL do not report useful feedback to help the transformation developers fix the fault. Consequently, manually examining the full MT and its contracts and reasoning on the implicit rule interactions become a time-consuming routine to debug MTs.

Because of the advancement in computer science in the last couple of decades (e.g. in the performance of satisfiability modulo theory -SMT -solvers), many researchers are interested in developing techniques that can partially or fully automate the localization of faults in software (we refer the reader to [START_REF] Roychoudhury | Formula-based software debugging[END_REF][START_REF] Wong | A survey on software fault localization[END_REF] for an overview).

In this work, we argue that the characteristics of the considered programming language have a significant impact on the precision and automation of the fault localization. More precisely, we think that in MT languages like ATL, automated fault localization can be more precise because of the available static trace information, i.e. inferred information among types of generated target elements and the rules that potentially generate these types. This idea has recently been introduced in [START_REF] Burgueño | Static fault localization in model transformations[END_REF] using a conservative and syntactical approach. However, we believe that a deductive approach can fully exploit its potential.

Our deductive approach is based on a set of sound natural deduction rules. It includes 4 rules for the ATL language based on the concept of static trace information, and 16 ordinary natural deduction rules for propositional and predicate logic [START_REF] Huth | Logic in Computer Science: Modelling and Reasoning About Systems[END_REF]. Then, we propose an automated proof strategy that applies these deduction rules on the input OCL postcondition to generate sub-goals. Each sub-goal contains a list of hypotheses deduced from the input postcondition, and a sub-case of the input postcondition to be verified. Successfully proving the sub-goals soundly implies the satisfaction of the input OCL postcondition. When a sub-goal is not verified, we exploit its hypotheses in two ways to help the user pinpoint the fault: (a) slicing the ATL MT into a simpler transformation context; (b) providing debugging clues, deduced from the input postcondition to alleviate the cognitive load for dealing with unverified sub-cases. Our fault localization approach has been implemented and integrated with VeriATL. We evaluate our approach with mutation analysis. The result shows that: (a) the guilty constructs are presented in the slice; (b) deduced clues assist developers in various debugging tasks (e.g. the elaboration of a counter-example); (c) the number of sub-goals that need to be examined to pinpoint a fault is usually small. Paper organization. We motivate our work by a sample problem in Section 2. Section 3 illustrates our fault localization approach in detail. Evaluation is presented in Section 4, followed by discussion of the limitations identified in our approach. Section 5 compares our work with related research, and Section 6 draws conclusions and lines for future work.

Motivating Example

As our running example we use the HSM2FSM MT. HSM2FSM transforms hierarchical state machine (HSM ) models to flattened state machine (FSM ) models. Both models conform to the same metamodel (Fig. 1). However, classifiers in the two metamodels are distinguished by the HSM and FSM prefix. Specifically, a named StateMachine contains a set of labelled Transitions and named AbstractStates. Each AbstractState has a concrete type, which is either RegularState, InitialState or CompositeState. A Transition links a source to a target AbstractState. Moreover, CompositeStates are only allowed in the models of HSM, and optionally contain a set of AbstractStates. Listing 1.1. The OCL contracts for HSM and FSM

Developing the ATL Transformation

Then, the developer implements the ATL transformation HSM2FSM (a snippet is shown in Listing 1.21 ). The transformation is defined via a list of ATL matched rules in a mapping style. The first rule maps each StateMachine element to the output model (SM2SM ). Then, we have two rules to transform AbstractStates: regular states are preserved (RS2RS ), initial states are transformed into regular states when they are within a composite state (IS2RS ). Notice here that initial states are deliberately transformed partially to demonstrate our problem, i.e. we miss a rule that specifies how to transform initial states when they are not within a composite state. The remaining three rules are responsible for mapping the Transitions of the input state machine.

Each ATL matched rule has a from section where the source pattern to be matched in the source model is specified. An optional OCL constraint may be added as the guard, and a rule is applicable only if the guard evaluates to true on the source pattern. Each rule also has a to section which specifies the elements to be created in the target model. The rule initializes the attributes/associations of a generated target element via the binding operator (<-). An important feature of ATL is the use of an implicit resolution algorithm during the target property initialization. Here we illustrate the algorithm by an example: 1) considering the binding stateMachine <-rs1.stateMachine in the RS2RS rule (line 7 of Listing 1.2), its right-hand side is evaluated to be a source element of type HSM!StateMachine; 2) the resolution algorithm then resolves such source element to its corresponding target element of type FSM!StateMachine (generated by the SM2SM rule); 3) the resolved result is assigned to the left-hand side of the binding. While not strictly needed for understanding this paper, we refer the reader to [START_REF] Jouault | ATL: A model transformation tool[END_REF] for a full description of the ATL language.

Formally Verifying the ATL Transformation

The source and target EMF metamodels and OCL contracts combined with the developed ATL transformation form a Hoare triple which can be used to verify the correctness of the ATL transformation, i.e. MM, Pre, Exec Post. The Hoare triple semantically means that, assuming the axiomatic semantics of the involved EMF metamodels (MM ) and OCL In previous work, Cheng et al. have developed the VeriATL verification system that allows such Hoare triples to be soundly verified [START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF]. Specifically, the VeriATL system describes in Boogie what correctness means for the ATL language in terms of structural Hoare triples. Then, VeriATL delegates the task of interacting with Z3 for proving these Hoare triples to Boogie. The axiomatic semantics of EMF metamodels and the OCL language are encoded as Boogie libraries in VeriATL. These libraries can be reused in the verifier designs of MT languages other than ATL.

In our example, VeriATL successfully reports that the OCL postcondition Post1 is not verified by the MT in Listing 1.2. This means that the transformation does not guarantee that each Transition has at least one source in the output model. Without any capability of fault localization, the developer then needs to manually inspect the full transformation and contracts to understand that the transformation is incorrect because of the absence of an ATL rule to transform InitialStates that are not within a CompositeState.

Our Goal: Localizing the Fault

In our running example, our proposed fault localization approach presents the user with two problematic transformation scenarios. One of them is shown in Listing 1.3. The scenario consists of the input preconditions (abbreviated at line 1), a slice of the transformation (abbreviated at lines 3 -5), and a sub-goal derived from the input postcondition. The sub-goal contains a list of hypotheses (lines 7 -12) with a conclusion (line 13).

The scenario in Listing 1.3 contains the following information, that we believe to be valuable in identifying and fixing the fault:

-Transformation slicing. The only relevant rules for the fault captured by this problematic transformation scenario are RS2RS, IS2RS and T2TC (lines 3 -5). -Debugging clues. The error occurs when a transition t0 is generated by the rule T2TC

(lines 8 -10), and when the source state of the transition is not generated (line 11). In addition, the absence of the source for t0 is due to the fact that none of the RS2RS and IS2RS rules is invoked to generate it (line 12).

From this information, the user could find a counter-example in the source models that falsifies Post1 (shown in the top of Fig. 2): a transition t c between an initial state i c (which is not within a composite state) and a composite state c c , where c c composites another initial state i c . This counter-example matches the source pattern of the T2TC rule (as shown in the bottom of Fig. 2). However, when the T2TC rule tries to initialize the source of the generated transition t2 (line 23 in Listing 1.2), i c cannot be resolved because there is no rule to match it. In this case, i c (of type HSM!InitialState)) is directly used to initialize the source of t2 (t2.source is expected to be a sub-type of FSM!AbstractState). This causes an exception of type mismatch, thus falsifying Post1. The other problematic transformation scenario pinpoints the same fault, showing that Post1 is not verified by the MT also when t0 is generated by T2TA. In the next section, we describe how we automatically generate problematic transformation scenarios like the one shown in Listing 1.3.

Overview of Fault Localization for ATL by Natural Deduction and Program Slicing

The flowchart in Fig. 3 shows a bird's eye view of our approach to enable fault localization for VeriATL. The process takes the involved metamodels, all the OCL preconditions, the ATL transformation and one of the OCL postconditions as inputs. We require all inputs to be syntactically correct. If VeriATL successfully verifies the input ATL transformation, we directly report a confirmation message to indicate its correctness (w.r.t. the given postcondition) and the process ends. Otherwise, we generate a set of problematic transformation scenarios (as the one shown in Listing 1.3), and a proof tree to the transformation developer.

To generate problematic transformation scenarios, we first perform a systematic approach to generate sub-goals for the input OCL postcondition. Our approach is based on a set of sound natural deduction rules. The set contains 16 rules for propositional and predicate logic such as introduction/elimination rules for ∧ and ∨ [START_REF] Huth | Logic in Computer Science: Modelling and Reasoning About Systems[END_REF], but also 4 rules Then, we design an automated proof strategy that applies the designed natural deduction rules on the input OCL postcondition. Executing our proof strategy generates a proof tree. The non-leaf nodes are intermediate results of deduction rule applications. The leafs in the tree are the sub-goals to prove. Each sub-goal consists of a list of hypotheses and a conclusion to be verified. The aim of our automated proof strategy is to simplify the original postcondition as much as possible to obtain a set of sub-conclusions to prove. As a by-product, we also deduce new hypotheses from the input postcondition and the transformation as debugging clues.

Next, we use the trace information in the hypotheses of each sub-goal to slice the input MT into simpler transformation contexts. We then form a new Hoare triple for each sub-goal consisting of the semantics of metamodels, input OCL preconditions, sliced transformation context, its hypotheses and its conclusion.

We send these new Hoare triples to the VeriATL verification system to check. Notice that successfully proving these new Hoare triples implies the satisfaction of the input OCL postcondition. If any of these new Hoare triples is not verified by VeriATL, the input OCL preconditions, the corresponding sliced transformation context, hypotheses and conclusion of the Hoare triple are presented to the user as a problematic transformation scenario for fault localization. The Hoare triples that were automatically proved by VeriATL are pruned away, and are not presented to the transformation developer. This deductive verification step by VeriATL makes the whole process practical, since the user is presented with a limited number of meaningful scenarios.

Then, the transformation developer consults the generated problematic transformation scenarios and the proof tree to debug the ATL transformation. If modifications are made on the inputs to fix the bug, the generation of sub-goals needs to start over. The whole process keeps iterating until the input ATL transformation is correct w.r.t. the input OCL postcondition.

Natural Deduction Rules for ATL

Our approach relies on 20 natural deduction rules (7 introduction rules and 13 elimination rules). The 4 elimination rules (abbreviated by X e ) that specifically involve ATL are shown in Fig. 4. The other rules are common natural deduction rules for propositional and predicate logic [START_REF] Huth | Logic in Computer Science: Modelling and Reasoning About Systems[END_REF]. Regarding the notations in our natural deduction rules:

-Each rule has a list of hypotheses and a conclusion, separated by a line. We use standard notation for typing (:) and set operations. -Some special notations in the rules are T for a type, M M T for the target metamodel, R n for a rule n in the input ATL transformation, x.a for a navigation expression, and i for a fresh variable / model element. In addition, we introduce the following auxiliary functions: cl returns the classifier types of the given metamodel, trace returns the ATL rules that generate the input type (i.e. the static trace information) 2 , genBy(i,R) is a predicate to indicate that a model element i is generated by the rule R, unDef(i) abbreviates i.oclIsUndefined(), and All(T) abbreviates T.allInstances().

x.a : T T ∈ cl(M M T )

x.a ∈ All(T ) ∨ unDef (x.a)

Tpe1

x.a :

Seq T T ∈ cl(M M T ) (|x.a| > 0 ∧ ∀i • (i ∈ x.a ⇒ i ∈ All(T ) ∨ unDef (i))) ∨ |x.a| = 0 Tpe2 T ∈ cl(M M T ) trace(T ) = {R 1 , ..., Rn } i ∈ All(T ) genBy(i, R 1 ) ∨ ... ∨ genBy(i, Rn) Tre1 T ∈ cl(M M T ) trace(T ) = {R 1 , ..., Rn } i : T unDef (i) ¬(genBy(i, R 1 ) ∨ ... ∨ genBy(i, Rn))
Tre2 Fig. 4. Natural deduction rules that specific to ATL Some explanation is in order for the natural deduction rules that are specific to ATL:

-First, we have two type elimination rules (T P e1 , T P e2 ). T P e1 states that every singlevalued navigation expression of the type T in the target metamodel is either a member of all generated instances of type T or undefined. T P e2 states that the cardinality of every multi-valued navigation expression of the type T in the target metamodel is either greater to zero (and every element i in the multi-valued navigation expression is either a member of all generated instances of type T or undefined) or equal to zero. -Second, we have 2 elimination rules for trace (T R e1 , T R e2 ). These rules state that, given that the rules R 1 ,...,R n in the input ATL transformation are responsible to create model elements of type T in the target metamodel, we may rightfully conclude that:

• (T R e1 ): every created element i of type T is generated by one of the rules R 1 ,...,R n .

• (T R e2 ): every undefined element i of type T is not generated by any of the rules R 1 ,...,R n .

Soundness of natural deduction rules. The soundness of our natural deduction rules is based on the operational semantics of the ATL language. Specifically, the soundness for type elimination rules T P e1 and T P e2 is straightforward. We prove their soundness by enumerating the possible states of initialized navigation expressions for target elements. Specifically, assuming that the state of a navigation expression x.a is initialized in the form x.a<-exp where x.a is of a non-primitive type T :

-If exp is not a collection type and cannot be resolved (i.e. exp cannot match the source pattern of any ATL rules), then x.a is undefined3 . -If exp is not a collection type and can be resolved, then the generated target element of the ATL rule that matches exp is assigned to x.a. Consequently, x.a could be either a member of All(T) (when the resolution result is of type T ) or undefined (when it is not). -If exp is of collection type, then all of the elements in exp are resolved individually, and the resolved results are put together into a pre-allocated collection col, and col is assigned to x.a.

The first two cases explain the two possible states of every single-valued navigation expressions (T P e1 ). The third case explains the two possible states of every multi-valued navigation expressions (T P e2 ). The soundness of trace elimination rules T R e1 is based on the surjectivity between each ATL rule and the type of its created target elements [START_REF] Büttner | On verifying ATL transformations using 'off-the-shelf' SMT solvers[END_REF]: elements in the target metamodel exist if they have been created by an ATL rule since standard ATL transformations are always executed on an initially empty target model. When a type can be generated by executing more than one rule, then a disjunction considering all these possibilities is made for every generated elements of this type.

About the soundness of the T R e2 rule, we observe that if a target element of type T is undefined, then clearly it does not belong to All(T). In addition, the operational semantics for the ATL language specifies that if a rule R is specified to generate elements of type T, then every target elements of type T generated by that rule belong to All(T) (i.e. R ∈ trace(T ) ⇒ ∀i • (genBy(i, R) ⇒ i ∈ All(T ))) [START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF]. Thus, T R e2 is sound as a logical consequence of the operational semantics for the ATL language (i.e. R ∈ trace(T ) ⇒ ∀i • (i / ∈ All(T ) ⇒ ¬genBy(i, R))).

Automated Proof Strategy

A proof strategy is a sequence of proof steps. Each step defines the consequences of applying a natural deduction rule on a proof tree. A proof tree consists of a set of nodes. Each node is constructed by a set of OCL expressions as hypotheses, an OCL expression as the conclusion, and another node as its parents node.

Next, we illustrate a proof strategy (Algorithm 1) that automatically applies our natural deduction rules on the input OCL postcondition. The goal is to automate the derivation of information from the postcondition as hypotheses, and simplify the postcondition as much as possible.

Algorithm 1 An automated proof strategy for VeriATL for each node leaf ∈ getLeafs(Tree) do 5:

Tree ← intro(leaf) ∪ Tree 6:

end for 7: while leafs = size(getLeafs(Tree)) 8: do 9:

leafs ← size(getLeafs(Tree)) 10:

for each node leaf ∈ getLeafs(Tree) do 11:

Tree ← elimin(leaf) ∪ Tree 12:

end for 13: while leafs = size(getLeafs(Tree))

Our proof strategy takes one argument which is one of the input postconditions. Then, it initializes the proof tree by constructing a new root node of the input postcondition as conclusion and no hypotheses and no parent node (line 1). Next, our proof strategy takes two sequences of proof steps. The first sequence applies the introduction rules on the leaf nodes of the proof tree to generate new leafs (lines 2 -7). It terminates when no new leafs are yield (line 7). The second sequence of steps applies the elimination rules on the leaf nodes of the proof tree (lines 8 -13). We only apply type elimination rules on a leaf when: (a) a free variable is in its hypotheses, and (b) a navigation expression of the free variable is referred by its hypotheses. Furthermore, to ensure termination, we enforce that if applying a rule on a node does not yield new descendants (i.e. whose hypotheses or conclusion are different from their parent), then we do not attach new nodes to the proof tree.

Transformation Slicing

Executing our proof strategy generates a proof tree. The leafs in the tree are the subgoals to prove by VeriATL. Next, we use the rules referred by the genBy predicates in the hypotheses of each sub-goal to slice the input MT into a simpler transformation context. We then form a new Hoare triple for each sub-goal consisting of the axiomatic semantics of metamodels, input OCL preconditions, sliced transformation context (Exec sliced ), its hypotheses and its conclusion, i.e. MM, Pre, Exec sliced , Hypotheses Conclusion.

If any of these new Hoare triples is not verified by VeriATL, the input OCL preconditions, the corresponding sliced transformation context, hypotheses and conclusion of the Hoare triple are constructed as a problematic transformation scenario to report back to the user for fault localization (as shown in Listing 1.3).

Our transformation slicing is based on the independence among ATL rules [START_REF] Tisi | Parallel execution of ATL transformation rules[END_REF]: each ATL rule is exclusively responsible for the generation of its output elements. Hence, when a sub-goal specifies a condition that a set of target elements should satisfy, the rules that do not generate these elements have no effects on the sub-goal. These rules can hence be safely sliced away.

Evaluation

In this section, we evaluate the practical feasibility and performance of our fault localization approach for the ATL language. The section concludes with a discussion of the obtained results and lessons learnt.

Research questions

We formulate two research questions to evaluate our fault localization approach:

(RQ1) Can our approach correctly pinpoint the faults in the given MT? (RQ2) Can our approach efficiently pinpoint the faults in the given MT?

Evaluation Setup

Our evaluation uses the VeriATL verification system [START_REF] Cheng | A sound execution semantics for ATL via translation validation[END_REF], which is based on the Boogie verifier (version 2.2) and Z3 (version 4.3). The evaluation is performed on an Intel 3 GHz machine with 8 GB of memory running the Windows operating system. VeriATL encodes the axiomatic semantics of the ATL language (version 3.7). The automated proof strategy and its corresponding natural deduction rules are currently implemented in Java.

To answer our research questions, we use the HSM2FSM transformation as our case study, and apply mutation analysis [START_REF] Jia | An analysis and survey of the development of mutation testing[END_REF] to systematically inject faults. In particular, we specify 14 preconditions and 5 postconditions on the original HSM transformation from [START_REF] Büttner | On verifying ATL transformations using 'off-the-shelf' SMT solvers[END_REF]. Then, we inject faults by applying a list of mutation operators defined in [START_REF] Burgueño | Static fault localization in model transformations[END_REF] on the transformation. We apply mutations only to the transformation because we focus on contract-based development, where the contract guides the development of the transformation. Our mutants are proved against the specified postconditions, and we apply our fault localization approach in case of unverified postconditions. We kindly refer to our online repository for the complete artifacts used in our evaluation [START_REF]A deductive approach for fault localization in ATL model transformations[END_REF].

Evaluation Results

Table 1 summarizes the evaluation results for our fault localization approach on the chosen case study. The first column lists the identity of the mutants4 . The second and third columns record the unverified OCL postconditions and their corresponding verification time. The fourth, fifth, sixth and seventh columns record information of verifying sub-goals, i.e. the number of unverified sub-goals / total number of sub-goals (4th), average verification time of sub-goals (5th), the maximum verification time among sub-goals (6th), total verification of sub-goals (7th) respectively. The last column records whether the faulty lines (L f aulty , i.e. the lines that the mutation operators operated on) are presented in the problematic transformation scenarios (P T S) of unverified sub-goals. First, we confirm that there is no inconclusive verification results of the generated subgoals, i.e. if VeriATL reports that the verification result of a sub-goal is unverified, then it presents a fault in the transformation. Our confirmation is based on the manual inspection of each unverified sub-goal to see whether there is a counter-example to falsify the sub-goal. This supports the correctness of our fault localization approach. We find that the deduced hypotheses of the sub-goals are useful for the elaboration of a counter-example (e.g. when they imply that the fault is caused by missing code as the case in Listing 1.3).

Second, as we inject faults by mutation, identifying whether the faulty line is presented in the problematic transformation scenarios of unverified sub-goals is also a strong indication of the correctness of our approach. Shown by the last column, all cases satisfies the faulty lines inclusion criteria. 3 out 10 cases are special cases (dashed cells) where the faulty lines are deleted by the mutation operator (thus there are no faulty lines). In the case of MF6#2, there are no problematic transformation scenarios generated since all the sub-goals are verified. By inspection, we report that our approach improves the completeness of VeriATL. That is the postcondition (#2) is correct under MF6 but unable to be verified by VeriATL, whereas all its generated sub-goals are verified.

Third, shown by the fourth column, in 5 out of 10 cases, the developer is presented with at most one problematic transformation scenario to pinpoint the fault. This positively supports the efficiency of our approach. The other 5 cases produce more sub-goals to examine. However, we find that in these cases each unverified sub-goal gives an unique phenomenon of the fault, which we believe is valuable to fix the bug. We also report that in rare cases more than one sub-goal could point to the same phenomenon of the fault. This is because the hypotheses of these sub-goals contain a semantically equivalent set of genBy predicates. Although they are easy to identify, we would like to investigate how to systematically filter these cases out in the future.

Fourth, from the third and fifth columns, we can see that each of the sub-goals is faster to verify than its corresponding postcondition by a factor of about 2. This is because we sent a simpler task than the input postcondition to verify, e.g. because of our transformation slicing, the Hoare triple for each sub-goal encodes a simpler interaction of transformation rules compared to the Hoare triple for its corresponding postcondition. From the third and sixth columns, we can further report that all sub-goals are verified in less time than their corresponding postcondition.

Limitations

Language coverage. In this work we consider a core subset of the ATL and OCL languages: (a) We consider the declarative aspect of ATL (matched rules) in non-refining mode, many-to-one mappings of (possibly abstract) classifiers with the default resolution algorithm of ATL. Non-recursive ATL helpers can be easily supported by inlining, and many-to-many mappings can be supported by extending our Boogie code generator. We also plan to investigate other features of ATL (e.g. lazy rules) to make our approach more general. (b) We support first-order OCL contracts and we plan to study more complex contracts in future work. Completeness of proof strategy. We define the completeness of a proof strategy meaning that every elements of target types referred by each sub-goal and every rule that may generate them are correctly identified after applying the proof strategy. If not detected, an incomplete proof strategy could cause our transformation slicing to erroneously slice away the rules that the sub-goal might depend on. By manual inspection, we confirm the completeness of our proof strategy in our case study. However, our proof strategy is in generally incomplete because: (a) we might lack deduction rules to continue the derivation of the proof tree; (b) our current proof strategy lacks of a backtracking mechanism when it chooses an unsuitable deduction rule to apply. Our current solution is detecting incomplete cases and reporting them to the user. In practice we check whether every elements of target types referred by each sub-goal are accompanied by a genBy predicate (this indicates full derivation). In future, we plan to improve the completeness of our approach by adding other natural deduction rules for ATL and smarter automated proof strategies.

Completeness of verification. Although we confirmed that there are no inconclusive sub-goals in our evaluation, our approach could report inconclusive sub-goals in general due to the underlying SMT solver. We hope the simplicity offered by the sub-goals would facilitate the user in making the distinction between incorrect and inconclusive sub-goals. In addition, if an input postcondition is inconclusive, our approach can help users to eliminate verified sub-goals to find the source of its inconclusiveness. Threats to validity of evaluation. We take a popular assumption in the fault localization community that multiple faults perform independently [START_REF] Wong | A survey on software fault localization[END_REF]. Thus, such assumption allows us to evaluate our fault localization approach in a one-postcondition-at-a-time manner. However, we cannot guarantee that this is the case for realistic and industrial MTs. We think classifying contracts into related groups could improve these situations. Scalability. The main scalability issue of our approach is that a complex OCL postcondition (e.g. an OCL expression with deeply nested quantifiers) can potentially generate a big number of sub-goals and corresponding problematic transformation scenarios. Verifying and displaying all of them becomes impractical for transformation developers. Since sub-goals are meant to be manually examined by the user, a reasonable solution is allowing the user to specify a bound for the maximum number of unverified cases to generate. To improve scalability we are also investigating the possibility of verifying intermediate nodes in the proof tree, and stop applying deduction rules if they are verified. Usability. Currently, our approach relies on the experience of transformation developer to interpret the deduced debugging clues. We think combining debugging clues with model finders would further help in debugging MT, e.g. by automatically generating the counterexamples [START_REF] Cuadrado | Uncovering errors in ATL model transformations using static analysis and constraint solving[END_REF].

Related Work

There is a large body of work on the topic of ensuring MT correctness [START_REF] Ab | A survey of approaches for verifying model transformations[END_REF]. To our knowledge our proposal is the first applying natural deduction with program slicing to increase the precision of fault localization in MT.

Büttner et al. use Z3 to verify a declarative subset of the ATL and OCL contracts [START_REF] Büttner | On verifying ATL transformations using 'off-the-shelf' SMT solvers[END_REF]. Their result is novel for providing minimal axioms that can verify the given OCL contracts. To understand the root of the unverified contracts, they demonstrate the UML2Alloy tool that draws on the Alloy model finder to generate counter examples [START_REF] Büttner | Verification of ATL transformations using transformation models and model finders[END_REF]. However, their tool does not guarantee that the newly generated counter example gives additional information than the previous ones. Oakes et al. statically verify ATL MTs by symbolic execution using DSLTrans [START_REF] Oakes | Fully verifying transformation contracts for declarative ATL[END_REF]. This approach enumerates all the possible states of the ATL transformation. If a rule is the root of a fault, all the states that involve the rule are reported.

Sánchez Cuadrado et al. present a static approach to uncover various typing errors in ATL MTs [START_REF] Cuadrado | Uncovering errors in ATL model transformations using static analysis and constraint solving[END_REF], and use the USE constraint solver to compute an input model as a witness for each error. Compared to their work, we focus on contract errors, and provide the user with sliced MTs and modularized contracts to debug the incorrect MTs.

Researchers have proposed several techniques that can partially or fully automate the localization of faults in software [START_REF] Roychoudhury | Formula-based software debugging[END_REF][START_REF] Wong | A survey on software fault localization[END_REF]. Program slicing refers to detect a set of program statements which could affect the values of interest [START_REF] Tip | A survey of program slicing techniques[END_REF][START_REF] Weiser | Program slicing[END_REF], which is used for fault localization of general programming languages. Few works have adapted this idea to localize faults in MTs. Aranega et al. define a framework to record the runtime traces between rules and the target elements these rules generated [START_REF] Aranega | Traceability mechanism for error localization in model transformation[END_REF]. When a target element is generated with an unexpected value, the transformation slices generated from the runtime traces are used for fault localization. While Aranega et al. focus on dynamic slicing, our work focuses on static slicing which does not require test suites to exercise the transformation.

The most similar approach to ours is the work of Burgueño et al. on syntactically calculating the intersection constructs used by the rules and contracts [START_REF] Burgueño | Static fault localization in model transformations[END_REF]. W.r.t. their approach we aim at improving the localization precision by considering also semantic relations between rules and contracts. This allows us to produce smaller slices by semantically eliminating unrelated rules from each scenario. Moreover, we provide debugging clues to help the user better understand why the sliced transformation causing the fault. However, their work considers a larger set of ATL. We believe that the two approaches complement each other and integrating them is useful and necessary.

We implement our approach in Java. However, we believe that integrating our approach to interactive theorem provers/framework such as Coq [START_REF] Bertot | Interactive Theorem Proving and Program Development: Coq'Art The Calculus of Inductive Constructions[END_REF] and Rodin [START_REF] Abrial | Rodin: An open toolset for modelling and reasoning in Event-B[END_REF] could be beneficial (e.g. drawing on recursive inductive reasoning). One of the easiest ways is through the Why3 language [START_REF] Filliâtre | Why3 -where programs meet provers[END_REF], which targets multiple theorem provers as its back-ends.

Conclusion and Future Work

In summary, in this work we confronted the fault localization problem for deductive verification of MT. We developed an automated proof strategy to apply a set of designed natural deduction rules on the input OCL postcondition to generate sub-goals. Each unverified subgoal yields a sliced transformation context and debugging clues to help the transformation developer pinpoint the fault in the input MT. Our evaluation with mutation analysis positively supports the correctness and efficiency of our fault localization approach. The result showed that: (a) faulty constructs are presented in the sliced transformation, (b) deduced clues assist developers in various debugging tasks (e.g. generate counter-example), (c) the number of sub-goals that need to be examined to pinpoint a fault are usually small.

Our future work includes facing the limitations identified during the evaluation (Section 4.4). We also plan to investigate how our decomposition can help us in reusing proof efforts. Specifically, due to requirements evolution, the MT and contracts are under unpredictable changes during the development. These changes can invalidate all of the previous proof efforts and cause long proofs to be recomputed. We think that our decomposition of sub-goals would increase the chances of reusing verification results, i.e. sub-goals that are not affected by the changes.

Fig. 1 .

 1 Fig. 1. The hierarchical and flattened state machine metamodel

Listing 1 . 3 .

 13 The problematic transformation scenario of the HSM2FSM transformation w.r.t. Post1 preconditions (Pre), by executing the developed ATL transformation (Exec), the specified OCL postcondition has to hold (Post).

Fig. 2 .

 2 Fig. 2. Counter-example derived from Listing 1.3 that falsify Post1

Fig. 3 .

 3 Fig. 3. Overview of providing fault localization for VeriATL

1 :

 1 Tree ← {createNode({}, Post, null)} 2: do 3: leafs ← size(getLeafs(Tree)) 4:

Table 1 .

 1 Evaluation metrics for the HSM2FSM case study

		Unveri. Post. ID Veri. Time(ms) Unveri. / Total Avg. Time (ms) Max Time (ms) Total Time (ms) Sub-goals	L f aulty ∈ P T S
	MT2 #5	3116	3 / 4	1616	1644	6464	True
	DB1 #5	2934	1 / 1	1546	1546	1546	-
	MB6 #4	3239	1 / 12	1764	2550	21168	True
	AF2 #4	3409	2 / 12	1793	2552	21516	True
	MF6	#2 #4	3779 3790	0 / 6 1 / 12	1777 1774	2093 2549	10662 21288	N/A True
	DR1	#1 #2	2161 2230	3 / 6 3 / 6	1547 1642	1589 1780	9282 9852	--
	AR	#1 #3	3890 4057	1 / 8 6 / 16	1612 1769	1812 1920	12896 28304	True True

Our HSM2FSM transformation is adapted from[START_REF] Büttner | On verifying ATL transformations using 'off-the-shelf' SMT solvers[END_REF]. The full version can be accessed at: https://goo. gl/MbwiJC.

In practice, we fill in the trace function by examining the output element types of each ATL rule, i.e. the to section of each rule.

In fact, the value of exp is assigned to x.a because of resolution failure. This causes a type mismatch exception and results in the value of x.a becoming undefined (we consider ATL transformations in nonrefinement mode where the source and target metamodels are different).

The naming convention for mutants are mutation operator Add(A) / Del(D) / Modify(M), followed by the mutation operand Rule(R) / Filter(F) / TargetElement(T) / Binding(B), followed by the position of the operand in the original transformation setting. For example, MB1 stands for the mutant which modifies the binding in the first rule.


S Tupin

J Molimard

V Cenizo

T Hoc

B Sohm

H Zahouani
email: hassan.zahouani@ec-lyon.fr

Multiscale approach to characterize mechanical properties of tissue engineered skin

Keywords: Reconstructed skin, Bi-photonic confocal microscopy, Indentation, Digital image correlation, Opto-mechanical system, Poisson's ratio

Tissue engineered skin usually consist of a multi-layered visco-elastic material composed of a fibrillar matrix and cells. The complete mechanical characterization of these tissues has not yet been accomplished. The purpose of this study was to develop a multiscale approach to perform this characterization in order to link the development process of a cultured skin to the mechanical properties. As a proof-of-concept, tissue engineered skin samples were characterized at different stages of manufacturing (acellular matrix, reconstructed dermis and reconstructed skin) for two different aging models (using cells from an 18-and a 61-year-old man). To assess structural variations, bi-photonic confocal microscopy was used. To characterize mechanical properties at a macroscopic scale, a light-load micro-mechanical device that performs indentation and relaxation tests was designed. Finally, images of the internal network of the samples under stretching were acquired by combining confocal microscopy with a tensile device. Mechanical properties at microscopic scale were assessed. Results revealed that adding cells during manufacturing induced structural changes, which provided higher elastic modulus and viscosity. Moreover, senescence models exhibited lower elastic modulus and viscosity. This multiscale approach was efficient to characterize and compare skin equivalent samples and permitted the first experimental assessment of the Poisson's ratio for such tissues.

Introduction

Skin tissue engineering has become an important field of research for medical and cosmetic applications. Most reconstructed skin models are validated by histological and biochemical characterization. The main limitation of existing models is their mechanical properties. Few mechanical studies have been performed that analyze the mechanical behavior of skin equivalents.

In vitro skin models are usually dermis-epidermis composites, which are based on threedimensional scaffolds that are colonized with dermal cells (fibroblasts) and covered with epidermal cells (keratinocytes). The 3D scaffold consists of either a natural polymer, synthetic polymers, or a mixture of both polymer types and is assembled to mimic the fibrillar structure of the dermis, which supports and anchors dermal cells. [START_REF] Cheung | A critical review on polymer-based bioengineered materials for scaffold development[END_REF][START_REF] Drury | Hydrogels for tissue engineering: scaffold design variables and applications[END_REF][START_REF] Mano | Natural origin biodegradable systems in tissue engineering and regenerative medicine: present status and some moving trends[END_REF] Therefore, the 3D scaffold is used as a template to guide cell growth and tissue development. Several parameters can be modulated to influence tissue formation, such as the polymer and its chemical composition, porosity and pore interconnectivity, 3D organization and mechanical properties. [START_REF] Kanungo | Density-property relationships in collagenglycosaminoglycan scaffolds[END_REF][START_REF] Peppas | Hydrogels in biology and medicine: from molecular principles to bionanotechnology[END_REF][START_REF] Shahabeddin | Characterization of skin reconstructed on a chitosan-cross-linked collagen-glycosaminoglycan matrix[END_REF] Because of the complex structure of skin, the mechanical behavior is mainly attributed to the structure of the dermal layer. [START_REF] Echinard | Under the epidermis, the dermis, or how to interpret skin cultures?[END_REF] Indeed, the fibrillar network, which is composed of collagen and elastin fibers, is linked together by physiological cross-linking and contributes to the resistance and elasticity of the tissue. [START_REF] Silver | Viscoelastic properties of human skin and processed dermis[END_REF][START_REF] Silver | Mechanobiology of force transduction in dermal tissue[END_REF] The surrounding medium, composed of proteoglycans and water, allows the fibers to slide against each other and therefore imparts viscous properties to the skin. The upper layer of the skin, the epidermis, has a stratified and cohesive structure resulting from the differentiation of keratinocytes. Engineered in vitro models have been used to demonstrate that the epidermal layer contributes to the mechanical properties of the entire skin when tensile tests are performed. [START_REF] Ebersole | Epidermal differentiation governs engineered skin biomechanics[END_REF] Skin equivalent tissues present a visco-elastic mechanical behavior. Many techniques, inspired from material science, have been developed to characterize the behavior of soft biological tissues, such as stress relaxation, [START_REF] Jachowicz | Indentometric analysis of in vivo skin and comparison with artificial skin models[END_REF] static indentation, [START_REF] Pailler-Mattei | In vivo skin biophysical behaviour and surface topography as a function of ageing[END_REF] dynamic indentation [START_REF] Boyer | Dynamic indentation on human skin in vivo: ageing effects[END_REF] and oscillating shear tests. [START_REF] Lamers | Large amplitude oscillatory shear properties of human skin[END_REF] Moreover, to characterize the viscous behavior of such biological materials, relaxation tests are often used. [START_REF] Hidaka | Elasticity and viscoelasticity of embolization microspheres[END_REF][START_REF] Pailler-Mattei | Rheological behaviour of reconstructed skin[END_REF] The optical visualization of the micro-structure is essential to understand the complexity of these biological samples. Only a few experiments have been performed that combined optics and mechanical testing. [START_REF] Boyer | Assessment of the in-plane biomechanical properties of human skin using a finite element model updating approach combined with an optical full-field measurement on a new tensile device[END_REF][START_REF] Gerhardt | A novel method for visualising and quantifying through-plane skin layer deformations[END_REF][START_REF] Roeder | Tensile mechanical properties of three-dimensional type I collagen extracellular matrices with varied microstructure[END_REF][START_REF] Voytik-Harbin | Simultaneous mechanical loading and confocal reflection microscopy for three-dimensional microbiomechanical analysis of biomaterials and tissue constructs[END_REF] In the present study, an original multiscale approach is developed to mechanically characterize tissue engineered skin and link its development process to the mechanical properties. Skin equivalent samples were characterized at different stages of manufacturing (acellular matrix, reconstructed dermis, reconstructed dermo-epidermis composite) for two different aging models (using cells from an 18-and a 61-year-old man). Fives samples of each type were produced. Bi-photonic confocal microscopy was firstly used to visualize the internal fibrillar and cellular network of the artificial samples. To characterize the global mechanical behavior of in vitro samples, a specific micro indentation/relaxation device was designed. Indentations were performed to characterize the elastic behavior of samples. The slope from the maximum of the unloading curve was used to determine the contact stiffness (𝐾 Z = d𝐹/d𝛿) and reduced Young's modulus (𝐸 * ) by following Hertz theory or the Johnson-Kendall-Roberts approach in the case of adhesion. Relaxation tests were performed on each sample. Experimental relaxation curves were analyzed using the generalized Maxwell model to characterize the instantaneous elastic modulus and the equivalent viscosity. Five indentation and five relaxation tests were performed on each sample (50 tests per skin equivalent type in total).

Finally, an opto-mechanical experimentation was developed by coupling bi-photonic confocal microscopy with a tensile device. It allowed us to acquire images of deformations in the micro-structure of the in vitro samples under unloaded and loaded (applied strain level of 2.5%) conditions. Digital image correlation was used to obtain displacement and strain fields for this specific area. It allowed the experimental assessment of the Poisson's ratio. After combining these data with the mechanical properties obtained by indentation, we were able to assess the internal stress and force on the fibers during tensile testing.

Materials and Methods

Material

Sample Types

This study was performed using reconstructed skin samples licensed by BASF ® Beauty Care Solutions, France. To evaluate the effect of sample structure on the mechanical properties, we mechanically characterized the acellular matrix, the reconstructed dermis and the reconstructed dermo-epidermis composite (reconstructed skin). Moreover, aging is a major concern in skin study. It is related to patient way of life, exposure to sun or pollution, and to his/her age. Here, in the context of reconstructed skin, it has been chosen to mimic the effect of aging by performing twice the entire process using cells from donors of different ages: an 18-and 61-year-old man. Fives samples of each type were produced.

Sample Preparation

The preparation of the reconstructed skin follows the work of Duplan-Perrat et al. [START_REF] Duplan-Perrat | Keratinocytes influence the maturation and organization of the elastin network in a skin equivalent[END_REF] and is summarized in Fig. 1.

Human keratinocytes and fibroblasts were isolated from human skin obtained from plastic surgery. Fibroblasts were grown in Dulbecco's modified Eagle's medium (DMEM with Glutamax-1, Life Technologies, Cergy-Pontoise, France) supplemented with 10% calf serum (HyClone, Logan, USA) and 100 µg/ml antibiotics (Normocin, InvivoGen, San Diego, USA). Keratinocytes were grown to confluence in defined K-SFM medium (Life Technologies) and antibiotics. Dermal equivalent fibroblasts were seeded at a density of 250,000 cells/cm 2 onto a dermal substrate consisting of a collagen-chitosan-glycosaminoglycan matrix (MIMEDISK™ BASF BCS, France). This reconstructed dermis (MIMEDERM™ BASF BCS, France) was grown for 28 days at 37 °C in a 5% CO2 atmosphere. The fibroblast DMEM medium (Life technologies) was supplemented with 10% fetal calf serum (FCS -HyClone, Logan, USA), 10 ng/ml epidermal growth factor, 50 µg/ml L-ascorbic acid

TUPIN et al. (2016)
(Sigma, Saint Louis USA), and antibiotics (Invivogen). The medium was changed every day.

Keratinocytes were seeded on the dermal equivalent on day 28 at a density of cells/cm 2 . Briefly, the cell-seeded matrices were submerged in keratinocyte medium consisting of a 3:1 mixture of DMEM and Ham's F12 (Invitrogen) supplemented with 10% calf serum (HyClone), 10 ng/ml epidermal growth factor (Sigma), 0.12 IU/ml insulin (Lilly, Saint-Cloud, France), 0.4 µg/ml hydrocortisone, 2 ng/ml triiodo-thyronine, 5 µg/ml apotransferrin, 0.4 µg/ml isoprenaline, 24.3 µg/ml adenine (Sigma) and antibiotics for days. Next, this reconstructed skin (MIMESKIN™ BASF BCS, France) was elevated to the air-liquid interface and cultured for 14 more days in a simplified keratinocyte medium containing DMEM supplemented with 0.8% bovine serum albumin (Sigma), 0.12 IU/ml insulin, 0.4 µg/ml hydrocortisone, 50 µg/ml L-ascorbic acid and antibiotics. The medium was changed every 2 days. Beauty Care Solutions, France.

Conservation Conditions

Nine samples of each type were used in the preset work. All of the following tests were performed in an air-conditioned room with controlled temperature (22-24°C) and relative humidity (30-40%). Prior to testing, the samples were stored in an incubator at a constant temperature of 37°C and 5% CO2. The tissue samples were removed from their culture plate and deposited on a glass slide for testing. Preliminary studies showed that mechanical properties are preserved at least 30 minutes; this time lapse is more than enough for the mechanical tests performed here.

Confocal Microscopy

A bi-photonic confocal microscope (Nikon, A1R MP Plus ® ) was used to acquire 3D images of the internal network of collagen and cells in our samples. Two visualization channels with an excitation of 820 nm were used. The first channel (500 to 550 nm in blue) visualized the auto-fluorescence (AF) of the cells. A second channel (400 to 492 nm in green) allowed us to visualize collagen using second harmonic generation (SHG). Stack of images were captured with slice dimensions of 1024 × 1024 pixels at a resolution of 0.5 µm/pixel and 1 µm between each slice. Density of signals were calculated for each XY slice as the percentage of pixels illuminated by each signal in the slice.

Mechanical Experimentations

Device

An original light-load indentation/relaxation device was developed for this type of test (previously described by Pailler-Mattei et al. [START_REF] Pailler-Mattei | Rheological behaviour of reconstructed skin[END_REF] ). The maximum force attained is 100 mN with a resolution of 0.02 mN. This device offers a displacement with a resolution of 3 nm and indenting velocities from 1 to 100 µm/s. A smooth PTFE spherical indenter with a radius of curvature R = 1.59 mm was selected as a probe to facilitate theoretical modeling and avoid adhesion between the probe and the biological sample.

Assessment of the Indentation Curve and Determination of the Stiffness and Young's Modulus

The principle of an indentation test is to push a rigid indenter perpendicularly into the sample in a load-unload displacement cycle at constant speed while the normal load, 𝐹 N , vs. the displacement, 𝛿, was recorded (Fig. 2). Before each indentation test, the zero contact point of indentation was estimated by positioning the indenter tip immediately above the specimen surface and moving the tip downwards at a velocity of 5 µm/s until a change in load greater than 0.02 mN was observed. The indenter was then raised by 1 mm to sever contact and moved down again to 200 µm above the surface. Finally, the force was zeroed.

After the force was zeroed, a force-controlled indentation test was conducted with an indentation speed of 25 µm/s and a maximum load of 0.5 mN. This maximum load was chosen such that the classical semi-infinite assumption in indentation data extraction remained valid (indentation depth of less than 10% of the sample height). For each sample, five indentations were performed and the load-displacement (N-m) data were recorded at a sampling of 1 kHz using specific software (based on LabVIEW ® ).

Following Zahouani et al. [START_REF] Zahouani | Characterization of the mechanical properties of a dermal equivalent compared with human skin in vivo by indentation and static friction tests[END_REF] , the contact stiffness, 𝐾 Z [N/m] is the slope of the initial portion of the unloading curve (Fig. 2) that is determined by:

𝐾 Z = | d𝐹 d𝛿 | 𝐹 max (1)
where 𝐹 [N] is the normal load and 𝛿 [m] is the displacement of the indenter. Following the Hertz theory, an approximation of the contact radius, 𝑎 [m] , depending on the maximum displacement of the indenter, 𝛿 [m] , and the indenter radius, 𝑅 [m] , is determined using Eq. ( 2):

𝑎 = √𝑅𝛿 (2)
Alternatively, the contact radius can by determined from the maximum load, 𝐹 max [N], contact stiffness, 𝐾 Z [N/m], and indenter radius, 𝑅 [m], as follows:

𝑎 = √ 3𝑅𝐹 max 2𝐾 Z (3)
Then, we can determine the reduced Young's modulus, 𝐸 * [Pa], using Eq. ( 4):

𝐸 * = 𝐾 Z 2𝑎 (4)
Note: The reduced Young's modulus is equal to:

1 𝐸 * = 1-𝜈 1 2 𝐸 1 + 1-𝜈 2 2 𝐸 2 (5)
where 𝐸 1 and 𝐸 2 are the Young's modulus of the sample and the Young's modulus of the spherical indenter, respectively, and ν 1 and ν 2 are the Poisson's ratio of the sample and the Poisson's ratio of the indenter, respectively. Therefore, 𝐸 2 » 𝐸 1 and the reduced Young's modulus is:

1 𝐸 * = 1-𝜈 1 2 𝐸 1
.

When the probe adhered to the sample, a negative force was observed (Fig. 2) and the slope of the unloading curve changed. Johnson et al. [START_REF] Johnson | Surface energy and the contact of elastic solids[END_REF] developed a successful elastic model, designated as the JKR model, which helped to correct for this adhesive behavior during estimation of the elastic modulus.

The new approximation of the contact area, 𝑎 JKR [m], is:

𝑎 JKR = √ 3𝑅 4𝐸 * (𝐹 max + 3𝜋𝑅Δ𝛾 + √6𝜋𝑅Δ𝛾𝐹 max + (3𝜋𝑅Δ𝛾) 2 ) 3 ( 6
)
where Δ𝛾 is the work of adhesion.

The maximal adhesion load can be expressed as:

𝐹 adh = - 3 2 𝜋𝑅Δ𝛾 (7)
Finally, the reduced elastic modulus can be obtained using Eq. ( 4), Eq. ( 6) and Eq. ( 7):

𝐸 JKR * = √ 𝐾 Z 3 6𝑅(𝐹 max -2𝐹 adh -2√𝐹 adh (𝐹 max -𝐹 adh )) (8)

Assessment of the Relaxation Curve and Determination of the Instantaneous Young's Modulus and Equivalent Viscosity

A relaxation test is similar to an indentation test with a lag time between the loading and unloading parts of the curve. During this time, we followed the decrease in the normal load, which is the relaxation curve (Fig. 3).

We followed an identical pre-test procedure for indentation using the following parameters: indentation speed of 50 µm/s, maximum load of 0.5 mN and latency time of 100 s. For each sample, five relaxation tests were conducted and load-time (N-s) data were recorded at a sampling of 1 kHz using developed software (on LabVIEW ® ).

The generalized Maxwell rheological model with two branches was applied to fit the relaxation part of the curves. Thanks to this model, the instantaneous elastic modulus, 𝐸 ins [Pa], and the equivalent viscosity, 𝜂 equ [Pa.s], can be calculated according to the contact area.

The generalized Maxwell model consists of a spring that is connected in parallel to nbranches that contain a spring and a dashpot connected in series (Fig. 3). Based on the generalized Maxwell model, an analytical solution using a functional equation for the spherical indentation load-relaxation test was proposed by Pailler-Mattei et al. [START_REF] Pailler-Mattei | Rheological behaviour of reconstructed skin[END_REF] :

𝐹(t) = 16 3 𝑅 1/2 𝛿 0 3/2 (𝐺 0 + ∑ 𝐺 𝑖 2 𝑖=1
e -𝑡/𝜏 𝑖 )

where 𝜏 𝑖 = 𝜂 𝑖 /𝐺 𝑖 . The experimental data obtained from the relaxation curves were fitted to the analytical models (Fig. 3) using the simplex method [START_REF] Nelder | A simplex algorithm for function minimization[END_REF] to obtain 𝐺 𝑖 and 𝜂 𝑖 . The Young's moduli, 𝐸 𝑖 , can be deduced directly from 𝐺 𝑖 :

𝐸 𝑖 = 2𝐺 𝑖 (1 + 𝜈 1 ) (10)
Finally, the instantaneous Young's modulus, 𝐸 ins , and equivalent viscosity, 𝜂 equ , can be obtained by reducing the model to a single branch:

𝐸 ins = 𝐸 0 + ∑ 𝐸 𝑖 2 𝑖=1
(11)

𝜂 eq = (∑ 𝐸 𝑖 𝑛 𝑖=1 ) 2 ∑ 𝐸 𝑖 𝑛 𝑖=1 /𝜏 𝑖 (12)

Opto-Mechanical Experimentation

Device

The goal of this experiment was to combine mechanical stress with an optical imaging technique to observe the response of the internal fibers to a macroscopic stress.

Combining the confocal microscope with a micro-tensile device allowed us to visualize the internal network of collagen under stretching. The behavior of the fibers during tensile testing was analyzed using digital image correlation.

The tensile device (DEBEN, 500N Compact Microtest ® ) was placed in situ in the confocal microscope. The following tensile parameters were used: a velocity of 2 mm/s and a maximum engineering strain of 2.5%.

We captured stacks of images with slice dimensions of 1024×1024 pixels at a resolution of 0.5 µm/pixel and 1 µm between each slice. The top of the dermis was analyzed, which is where collagen is synthesized. Each stack (∼60 µm in Z) was transformed into a single image using maximum-intensity projection (ImageJ ® software). This image represents the condition of the area at a specific macroscopic strain. Finally, we obtained two images per sample that displayed internal views of the samples at their initial state and after a strain of 2.5%.

Digital Image Correlation

The developed method is based on digital image correlation (DIC) [START_REF] Sutton | Image Correlation for Shape, Motion and Deformation Measurements-Basic Concepts, Theory and Applications[END_REF] and was previously described by Molimard et al. [START_REF] Molimard | Frequency-based image analysis of random patterns: an alternative way to classical stereocorrelation[END_REF] Only the basic principle is given here.

2D Digital Image Correlation

Assuming a reference image 𝑖𝑚 0 described by 𝑓(𝑟, 𝑠), a deformed image 𝑖𝑚 1 of 𝑖𝑚 0 after a small deformation is described by 𝑔(𝑟, 𝑠):

𝑔(𝑟, 𝑠) = 𝑓(𝑟 -𝛿 𝑥 , 𝑠 -𝛿 𝑦 ) + 𝑏(𝑟, 𝑠)

where 𝛿 𝑥 and 𝛿 𝑦 are the components of the displacement of 𝑖𝑚 1 and 𝑏(𝑟, 𝑠) is the noise. One method to determine 𝛿 𝑥 and 𝛿 𝑦 is to maximize function ℎ, which is defined by:

ℎ(𝑟, 𝑠) = (𝑔 ⊗ 𝑓)(𝑟, 𝑠) = ∫ ∫ 𝑔(𝑎, 𝑏)𝑓(𝑎 -𝑟, 𝑏 -𝑠) +∞ -∞ +∞ -∞ d𝑎 d𝑏 (14)
where ⊗ denotes the cross-correlation product. The position of the maximum value of ℎ(𝑟, 𝑠) corresponds to the maximum displacement probability (𝛿 𝑥 , 𝛿 𝑦 ).

Equation 14 is used at a local scale on the region of interest (ROI). To provide a final displacement field, the procedure is repeated for regions located all over the map.

Sub-pixel Algorithm

The sub-pixel algorithm is based on a phase analysis of each ROI before and after mechanical stimulation. Each phase is proportional to the displacement normal to the corresponding direction. In the absence of a phase jump, the displacements can be related to any defined phase using the following relationship:

{ ⋮ Δ𝜙 𝜃 𝑖 ⋮ } = [ ⋮ ⋮ 2𝜋 𝑝 𝜃 𝑖 𝑐𝑜𝑠𝜃 𝑖 2𝜋 𝑝 𝜃 𝑖 𝑠𝑖𝑛𝜃 𝑖 ⋮ ⋮ ] = 𝐴 { 𝛿 𝑥 𝛿 𝑦 } (15)
where{ ⋮ Δ𝜙 𝜃 𝑖 ⋮ } is a vector containing all of the phase differences, { 𝛿 𝑥 𝛿 𝑦 } is the displacement vector, and 𝐴 is the sensitivity matrix. 𝐴 is a rectangular matrix with 2 rows and as many lines as defined phases, which can be up to 255 under usual conditions.

The displacements can then be derived from Eq.( 15) using the pseudo-inverse of 𝐴. This operation is possible if det((𝐴 𝑡 𝐴)) ≠ 0. In practice, this indicates that at least two phases along two different directions exist.

{ 𝛿 𝑥 𝛿 𝑦 } = (A t A) -1 A t { ⋮ Δ𝜙 𝜃 𝑖 ⋮ } (16)

Strain Fields

Finally, the strain fields ( 𝜀 𝑋𝑋 and 𝜀 𝑌𝑌 ) were obtained using a 2D local polynomial approach. This method was proven to be 8 times more efficient than classical finite difference methods. 1

Global Strain Values and Poisson's Ratio Assessment

The distribution of strain fields was supposed to follow a normal distribution, and the global strain values were evaluated as the mean. Using those values, we determined the Poisson's ratio, 𝜈, which is defined for a tensile test as:

𝜈 = -𝜀 𝑌𝑌 𝜀 𝑋𝑋 (17)

Determination of Mechanical Properties at the Microscopic Scale

Following linear elasticity and small deformation hypotheses, we can determine the mechanical properties at a local level. Here, considering the microscope magnification, the local level is a microscopic one, corresponding to the fiber scale.

Using the global strain values, the Poisson's ratio and the elastic modulus, which was characterized by indentation and corrected with the experimental value of the Poisson's ratio, we determined the stress along 𝑥 → , 𝜎 𝑋𝑋 [Pa], which is defined for a tensile test as:

𝜎 𝑋𝑋 = 𝐸 1-𝜈 2 (𝜀 𝑋𝑋 + 𝜈𝜀 𝑌𝑌 ) (18)
The work of elastic deformation along 𝑥 → , 𝑊 𝑒𝑋 [J], can be determined using the volume of optical acquisition, 𝑉 [m 3 ], as follows:

𝑊 𝑒𝑋 = 1 2 𝐸𝜀 𝑋𝑋 2 𝑉 (19)
Finally, using the first theorem of Castigliano (Eq.20), forces along 𝑥 → , 𝐹 𝑋 [N], that are applied to the fibers in this volume are determined by Eq. ( 21):

𝑊 𝑒𝑋 = 1 2 𝐹 𝑋 𝛿 𝑥 ( 20
)
𝐹 𝑋 = 2𝑊 𝑒𝑋 𝛿 𝑥 (21)

Statistics

Mean and standard deviation values were calculated for the results of the same tissue sample group. The results are presented as M ± SD with M = mean and SD = standard deviation. Moreover, unpaired Student's t tests were performed between two groups, and the p values, p, are presented.

Results

Micro-structure Observations

Acellular Matrix

The acellular matrix, basis of next samples, is shown in Fig. 8a. It is a collagen-chitosanglycosaminoglycan scaffold which could be visualized in both channels thanks to the second harmonic generation of collagen and the auto-fluorescence of chitosan, previously described in the literature. 4

Reconstructed Dermis

Confocal images of reconstructed dermis samples were acquired (Fig. 4). Because two channels were used, we were able to visualize the signal of the artificial matrix (visible by auto-fluorescence of the chitosan in blue) and the signal of the human collagen (visible by second harmonic generation of the collagen in green). From these images we observed that the fibroblasts induced the neo-synthesis of collagen. Moreover, the cells from the younger man synthesized more collagen. Indeed, the density of synthetized collagen (in green) is equivalent at the top of the sample but decrease strongly in depth with the use of the older man's cells.

Reconstructed Skin

Confocal images were acquired for the reconstructed skin sample (Fig. 5). At the bottom of the images, we can observe the matrix (visible by auto-fluorescence of the chitosan in blue) and collagen neo-synthesis (visible by second harmonic generation in green). Keratinocytes were seeded on top of this layer. As described in the literature, these cells increase the production of keratin during their differentiation. [START_REF] Wu | Depth-resolved fluorescence spectroscopy reveals layered structure of tissue[END_REF] This produced a highly fluorescent region that corresponded with the keratinizing layer (stratum corneum) on the surface of the sample. Keratinocytes in the lower section of this reconstructed epidermis produced less keratin; they are not visible in these images because their fluorescent intensity is significantly lower than that of the stratum corneum. The cells from the younger man synthesized more collagen. Indeed, the density of synthetized collagen (in green) is lower with the use of the older man's cells (30 and 18%, with the use of the younger man's cells and of the older man's cells, respectively).

Mechanical Results

Indentation

Reduced elastic modulus obtained from the indentation tests are presented in Fig. 6.

The results reveal increasing values (p < 0.001) from the acellular matrix (14.8 ± 1.1 kPa) to the reconstructed dermis (59.1 ± 6.6 kPa) and finally the reconstructed skin (160.3 ± 32.5 kPa). We were able to follow the effect of adding cells during the culture of artificial skin. Moreover, the use of cells from the older donor yielded significantly lower values of the reduced elastic modulus (-13.1 kPa, p < 0.01, and -76.6 kPa, p < 0.001, for the dermis and skin, respectively).

Relaxation

The mechanical parameters obtained from the relaxation tests are presented in Fig. 7.

The results exhibit a significant increase (p < 0.001) in the instantaneous elastic modulus from the acellular matrix (3.8 ± 0.3 kPa) to the reconstructed dermis (8.9 ± 1.3 kPa) and the reconstructed skin (83.3 ± 33.1 kPa). We are able to monitor the effect of adding cells while culturing the artificial skin. Moreover, the use of cells from the older donor yielded a lower instantaneous elastic modulus (-5.0 kPa, p < 0.001, and -33.3kPa, p < 0.01, for the dermis and skin, respectively).

In terms of viscosity, 𝜂 equ follows the same trend as the instantaneous elastic modulus. The viscosity increased from 2.5 ± 1.0 kPa.s to 7.1 ± 2.3 kPa.s and finally to 150.8 ± 97.5 kPa.s for the acellular matrix, reconstructed dermis, and reconstructed skin, respectively. A decrease in the equivalent viscosity was observed for the samples using older male cells (-2.5 kPa.s, p < 0.01, and -67.9 kPa.s, p < 0.05, for the dermis and skin, respectively).

Displacement Fields and Poisson's Ratio

For each sample, digital image correlation was applied to the top part of the matrix. An example of this image correlation process on the artificial matrix can be followed in Fig. 8.

Displacement fields (Fig. 8c and8d) are typical results for a tensile test. After derivation, strain maps (Fig. 8e and8f) were assessed and the global strains were estimated. Global strains are estimated by fitting the distribution of strain maps with a normal law. Using these parameters, the Poisson's ratio could be experimentally assessed (Table 1). The results exhibit a decrease in the Poisson's ratio from the acellular matrix (0.42) to the reconstructed dermis (0.21) and the reconstructed skin (0.20).

Moreover, the use of cells from the older donor did not yield a significant difference.

Mechanical Properties at the Microscopic Scale

Using previous parameters (global strains, Poisson's ratio and elastic modulus estimated by indentation), mechanical properties were calculated at the microscopic scale (Table .  2).

The results exhibit an increase in the internal stress and force from the acellular matrix (239.0 Pa and 11.6 µN, respectively) to the reconstructed dermis (948.6 Pa and 52.9 µN, respectively) and the reconstructed skin (2278 Pa and 130.4 µN, respectively). We are able to monitor the effect of adding cells while culturing the artificial skin.

Moreover, the use of cells from the older donor yielded lower parameters (-389.9 Pa and -25.6 µN for the reconstructed dermis and -1238 Pa and -74.3 µN for the reconstructed skin, respectively).

Structural and Aging Effects

Our work allowed us to study the effect of aging on artificial skin samples using cells from two donors of different ages (18-and 63-year-old). Mechanical experiments provided information on the elasticity and viscosity of those samples. Using a collagen-chitosanglycosaminoglycan matrix, we can create a reconstructed dermis using fibroblasts. This scaffold is composed of large fibers and is porous enough to be seeded with fibroblasts. After a period of culture, the cells colonized and began to remodel the matrix. The cells degraded the matrix and stimulated neo-synthesis of human collagen, as observed with bi-photonic confocal microscopy. With regards to the mechanical properties, the elastic modulus increased significantly between the matrix and the reconstructed dermis. This result was confirmed by two mechanical tests (indentation & relaxation). The use of older cells yielded a smaller elastic modulus that was equivalent to the initial matrix. This result is attributed to the less efficient collagen synthesis that was observed for older cells (Fig. 4b).

To reconstruct a layer of skin, keratinocytes were seeded on top of the reconstructed dermis. As observed with confocal microscopy, these cells formed a layer on top of the sample, which is the result of stratification and the complete differentiation of keratinocytes. The skin sample is composed of a rigid layer on top of a fibrillar milieu.

With regards to mechanical properties, the elastic modulus and viscosity increased significantly with the addition of cells. The rigid layer of differentiated keratinocytes maintained the entire structure of the sample, which increased the viscosity. Confocal images revealed similar layers of keratinocytes for matrices seeded with either young or old cells. Moreover, seeding keratinocytes on top of the reconstructed dermis stimulated fibroblast metabolism: the synthesis of collagen seemed to be equivalent for young and old cells (Fig. 5a and5b).

However, under the epidermis, the layer of fibroblasts was larger for matrices seeded with cells from the older donor. This structural configuration had a lower elastic modulus, as shown in Fig. 7.

Opto-mechanical Assessments

We know that in vivo skin aging induces a loss of collagen and elastic fiber density that decreases internal tension, as revealed by Pailler-Mattei et al. [START_REF] Pailler-Mattei | In vivo skin biophysical behaviour and surface topography as a function of ageing[END_REF] The goal of our optomechanical experiments was to reproduce this internal tension using in vitro samples.

The microstructure of the sample, before and after a macroscopic strain of 2.5%, was analyzed using digital image correlation to assess the displacement and strain fields (after 2D local polynomial derivation). It leads to the evaluation of the Poisson's ratio and the internal stress and force in the sample. This internal force reflects the ability of the sample to resist deformation.

The Poisson's ratio could be experimentally assessed (Table 1). The acellular matrix exhibits a value of 0.42 which is not far away from the usually used value for an incompressible material (0.45-0.5). The addition of cells in the acellular matrix induces a strong decrease of this parameter to a value of 0.2. This value tends to the experimental values assessed in vivo by Boyer et al. [START_REF] Boyer | Assessment of the in-plane biomechanical properties of human skin using a finite element model updating approach combined with an optical full-field measurement on a new tensile device[END_REF] Internal stress and force increased with the addition of cells during the manufacturing of the reconstructed skin (Table 2). This result emphasizes the participation of fibroblastsynthesized collagen during dermis culture and the addition of the stratified epidermis, which is composed of differentiated keratinocytes, in the manufacturing of reconstructed skin. Moreover, the results revealed a decrease in the internal stresses and forces when cells from the older male donor were used. This is attributed to the less efficient synthesis of collagen, which was visualized using bi-photonic confocal microscopy.

Conclusion

Measuring the visco-elastic properties and the Poisson's ratio of artificial skin samples are key initiatives to better understand their mechanical behavior.

Complete device and methods for the in vitro assessment of those properties at macroscopic scale were presented (indentation and relaxation tests). The results were obtained for artificial skin samples developed by BASF ® Beauty Care Solutions France; a variety of skin models were selected, based on manufacturing steps, and analyzed. The results clearly demonstrate the effect of those steps on the mechanical properties of the skin substitutes at macroscopic level. Moreover, the effect of aging was studied using senescence models (cells from the older male donor).

Figure 1 :

 1 Figure 1: Preparation of the tissue engineered samples provided by BASF ® Beauty Care Solutions, France.

Figure 2 :

 2 Figure 2: Experimental indentation curve of a sample presenting a visco-elastic adhesive behavior and schematic of the indentor penetrating the sample.

Figure 3 :

 3 Figure 3: Experimental relaxation curve of a sample presenting a visco-elastic behavior (in blue) and fitting curve using the generalized Maxwell model with two branches (in red).

Figure 4 :

 4 Figure 4: ZX view of reconstructed dermis samples obtained using bi-photonic confocal microscopy and plots of the density of signals in each image through depth.

Figure 5 :

 5 Figure 5: ZX view of reconstructed skin samples obtained using bi-photonic confocal microscopy and plots of the density of signals in each image through depth.

Figure 6 :

 6 Figure 6: Reduced elastic modulus, E*, obtained by indentation test on the tissue engineered samples. ( * p<0.05)

Figure 7 :

 7 Figure 7: Instantaneous elastic modulus, Eins, and equivalent viscosity, ηequ, obtained by relaxation test on the tissue engineered samples. ( * p<0.05)

Figure 8 :

 8 Figure 8: Digital image correlation applied to an acellular matrix sample: biphotonic confocal images (initial state (a) and after a strain of 2.5% (b)), displacement fields along X (c) and Y (d) evaluated using digital image correlation and strain fields εXX (e) and εYY (f) after derivation.

Table 1 : Global strains, εXX and εYY, and Poisson's coefficient, 𝝂, estimated for each tissue engineered sample.

 1

	Sample	𝜀 𝑋𝑋	𝜀 𝑌𝑌	𝜈
	Acellular matrix 0.016 -0.006 0.42
	Dermis 18yo	0.016 -0.003 0.21
	Dermis 61yo	0.012 -0.003 0.25
	Skin 18yo	0.014 -0.003 0.20
	Skin 61yo	0.012 -0.002 0.20

Table 2 : Mechanical parameters at microscopic level for each tissue engineered sample: stress, σX, elastic work, WeX, and force, FX.

 2 Sample𝜎 𝑋 [Pa] 𝑊 𝑒𝑋 [nJ] 𝐹 𝑋[µN]

	Acellular matrix 239.0	0.312	11.6
	Dermis 18yo	948.6	1.20	52.9
	Dermis 61yo	558.7	0.534	27.3
	Skin 18yo	2278	2.61	130.4
	Skin 61yo	1040	1.02	56.1
	4 Discussion

Confocal microscopy images were acquired to better understand the structural variations between samples. The images clearly show the effect of the cell type on collagen synthesis. Indeed, cells from the older man were less efficient at synthesizing collagen, which explains the decreased mechanical properties.

Finally, opto-mechanical experiments were performed, which combined bi-photonic confocal microscopy with a tensile device. The micro-structure, before and after 2.5% strain, was compared using digital image correlation. The displacement and strains that were obtained were used to estimate the Poisson's ratio of the samples. The values of this parameter for our skin equivalent samples vary between 0.2 and 0.25. Also, mechanical properties at the microscopic scale were evaluated. This experiment demonstrates how much a sample can resist a specific strain. Senescence models are less efficient than normal ones, which can be correlated to in vivo observations for the loss of tension in the skin that is observed with aging.

Acknowledgements

We would like to convey our deepest gratitude to BASF ® Beauty Care Solutions France (F-69007 Lyon) for providing us with the in vitro samples.

We would also like to thank the IVTV (ANR-10-EQPX-06-01) team for their help during imaging.

Conflict of interest

No conflicts of interest present in this study.