Frances M Davis

Yuanming Luo

Stéphane Avril

Ambroise Duprey

Jia Lu
email: jia-lu@uiowa.edu

Pointwise characterization of the elastic properties of planar soft tissues: Application to ascending thoracic aneurysms

Keywords: bulge inflation test, heterogeneous material properties, digital image correlation, inverse elastostatic analysis, thoracic aneurysm

In this manuscript, we present a combined experimental and computational technique that can identify the heterogeneous elastic properties of planar soft tissues. By combining inverse membrane analysis, digital image correlation, and bulge inflation tests, we are able to identify a tissue's mechanical properties locally. To show how the proposed method could be implemented, we quantified the heterogeneous material properties of a human ascending thoracic aortic aneurysm (ATAA). The ATAA was inflated at a constant rate using a bulge inflation device until it ruptured. Every 3 kPa images were taken using a stereo digital image correlation system. From the images, the three dimensional displacement of the sample surface was determined. A deforming NURBS mesh was derived from the displacement data and the local strains were computed. The wall stresses at each pressure increment were determined using inverse membrane analysis. The local material properties of the ATAA were then identified using the pointwise stress and strain data. To show that it is necessary to consider the heterogeneous distribution of the mechanical properties in the ATAA, three different forward finite element simulations using pointwise, elementwise, and homogeneous material properties were compared. The forward finite element predictions revealed that heterogeneous nature of the ATAA must be

Introduction

Most soft biological tissues are heterogeneous, nonlinear, and anisotropic. Yet, these heterogeneous tissues are often modeled as homogeneous. For some tissues the assumption of homogeneity is justified but, frequently the assumption of homogeneity is simply convenient. The identification of the heterogeneous material properties can be particularly important when attempting to understand the interaction between disease progression or tissue repair and a tissue's mechanical properties. In, for example, aneurysm rupture and growth, the identification of a tissue's heterogeneous properties is essential. In aortic aneurysms, spatial variations in the mechanical properties occur naturally due to local remodeling of the vessel wall [START_REF] Humphrey | Dysfunctional mechanosensing in aneurysms[END_REF][START_REF] Wilson | Importance of initial aortic properties on the evolving regional anisotropy, stiffness and wall thickness of human abdominal aortic aneurysms[END_REF]). This regional variation in the mechanical behavior has a direct effect on the wall stress, which is the most popular biomechanical factor for assessing aneurysm rupture risk [START_REF] Vorp | Biomechanics of abdominal aortic aneurysm[END_REF]. Furthermore, studying how the spatial distribution of the mechanical properties evolves as the aneurysm geometry and composition change could provide significant insight into the growth and remodeling process. Although the heterogeneous nature of aortic aneurysms is recognized, only a few investigators have characterized the spatial variation in mechanical properties [START_REF] Choudhury | Local mechanical and structural properties of healthy and diseased human ascending aorta tissue[END_REF][START_REF] Haskett | Microstructural and biomechanical alterations of the human aorta as a function of age and location[END_REF][START_REF] Iliopoulos | Regional and directional variations in the mechanical properties of ascending thoracic aortic aneurysms[END_REF][START_REF] Annaidh | Characterization of the anisotropic mechanical properties of excised human skin[END_REF][START_REF] Romo | In vitro analysis of localized aneurysm rupture[END_REF].

The primary barrier to investigating the heterogeneous nature of soft tissues is the difficulty developing experimental protocols that can accurately capture the spatial variation in mechanical properties. Conventional testing approaches, such as uniaxial and biaxial tensile tests, rely on the assumption of uniformity of the stress and strain in the center of the specimen. The most popular strategy for studying heterogeneity in soft tissues involves dividing the tissue into multiple test samples and characterizing the homogeneous response of each sample using uniaxial or biaxial tests. With this approach the regional heterogeneity of a tissue can be studied at a spatial resolution of centimeters. Using uniaxial or biaxial tests to study spatial variation in the material properties has two significant limitations: It requires the heterogeneity to be rather large sized and a significant number of samples.

The objective of this manuscript is to advance our ability to characterize heterogeneities in soft tissues in vitro. With this in mind, we restrict ourselves to planar soft tissues that can be modeled as membranes. From the mechanics perspective, membranes are thin structures that offer little resistance to bending. Many soft tissues can be modeled as membranes such as aneurysms, the bladder, skin, and the sclera in the eye [START_REF] Humphrey | Review Paper: Continuum biomechanics of soft biological tissues[END_REF][START_REF] Myers | The inflation response of the posterior bovine sclera[END_REF][START_REF] Tonge | Full-field bulge test for planar anisotropic tissues: Part I-Experimental methods applied to human skin tissue[END_REF]. Recently, there have been significant computational and theoretical developments in the formulation of inverse approaches to describe the finite deformation of membranes [START_REF] Lu | Pointwise identification of elastic properties in nonlinear hyperelastic membranes-Part I: Theoretical and computational developments[END_REF][START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF][START_REF] Zhao | Pointwise identification of elastic properties in nonlinear hyperelastic membranes -Part II: Experimental validation[END_REF]. To characterize the heterogeneous distribution of material properties in planar biological tissues we exploit these developments and the unique properties of membranes.

There are three primary components to the material characterization: bulge inflation tests, stereo digital image correlation (DIC), and inverse membrane analysis. By coupling bulge inflation tests with DIC to obtain full-field displacement measurements, we are able to compute the strain locally, as opposed to globally with traditional experimental protocols. The inverse membrane formulation capitalizes on a unique property of the membrane equilibrium equations, static determinacy, that allows the stress field to be identified without prior knowledge of the membrane's material properties [START_REF] Genovese | Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder in vitro[END_REF][START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF][START_REF] Romo | In vitro analysis of localized aneurysm rupture[END_REF]. A straight-forward pointwise characterization of a membrane's elastic properties can then be completed using the local stress-strain data from the inverse membrane analysis and the DIC measurements [START_REF] Genovese | Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder in vitro[END_REF][START_REF] Lu | Pointwise identification of elastic properties in nonlinear hyperelastic membranes-Part I: Theoretical and computational developments[END_REF][START_REF] Zhao | Pointwise identification of elastic properties in nonlinear hyperelastic membranes -Part II: Experimental validation[END_REF]Zhao et al., , 2011a,b),b). In this manuscript, we describe a protocol for the pointwise identification of anisotropic material properties in planar soft tissues. The method is applied to a human ascending thoracic aortic aneurysm (ATAA). The ATAA was subjected to a bulge inflation test and the displacement and pressure data collected are used as inputs to the inverse membrane method. The wall stresses were calculated using an inverse membrane formulation with a NURBS representation [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF]. Using the computed full field stress-strain data, the pointwise distribution of the mechanical properties for the ATAA specimens were identified. To substantiate the claim that a pointwise resolution is necessary to characterize the mechanical properties of the ATAA a forward finite element analysis was performed using homogeneous, elementwise, and pointwise distributions of the mechanical properties and the resulting errors in the predicted displacement, strain, and stress are reported.

Methods

Experimental Protocol

An unruptured ATAA section was collected from a patient undergoing elective surgery to replace his ATAA with a graft in accordance with a protocol approved by the Institutional Review Board of the University Hospi-tal Center of St. Etienne. After retrieval, the specimen was placed in saline solution and stored at 4 • C until testing, which occurred within 48 hours of the surgery. Immediately prior to testing, the ATAA was cut into a square specimen approximately 45 x 45 mm. Any fatty deposits were removed from the surface of the tissue to ensure that during mechanical testing the tissue did not slip in the clamps. An average thickness was found for the sample by measuring the thickness of the tissue at a minimum of 5 locations.

The specimen was clamped in the bulge inflation device, Fig. 1, so that the luminal side of the tissue faced outward. Then a speckle pattern was applied to the luminal surface using black spray paint. The sample was inflated using a syringe pump driven at 2 mL/min to infuse water into the cavity behind the sample. During the test, the pressure was measured using a digital manometer (WIKA, DG-10). Images of the inflating specimen were collected using a commercial DIC system (GOM, 5M LT) composed of two 8-bit CCD cameras equipped with 50 mm lenses (resolution: 1624 x 1236 px). The cameras were positioned 50 cm apart at an angle of 30 • with an aperture of f /11. This produced a depth of field of 15.4 mm which was sufficient to capture the deformation of the tissue up to failure. Images of the deforming sample were collected every 3 kPa until the sample ruptured.

After rupture, the collected images were analyzed using the commercial correlation software ARAMIS -(GOM, v. 6.2.0) to determine the three dimensional displacement of the tissue surface. For the image analysis, a facet size of 21 px and a facet step of 5 px were chosen based on the speckle pattern dot size, distribution, and contrast. The selected parameters produced a cloud of approximately 15,000 points where the three displacement values were calculated.

A preliminary study was performed to estimate the displacement and strain noise floor for the DIC measurements. Four samples were inflated to 0, 3, 15, 30, 45, and 60 kPa and a baseline stage was captured and 1 s later another stage was captured under the same loading conditions. The baseline image was correlated with the second stage to determine the displacement at every point i. The error, l, was calculated at each image point [START_REF] Myers | The inflation response of the posterior bovine sclera[END_REF]

l i = (r x i -R x i ) 2 + (r y i -R y i ) 2 + (r z i -R z i ) 2 (1)
where r x i , r y i , and r z i are the components of displacement in the second image and R x i , R y i , and R z i are the components of displacement in the baseline image. The mean displacement error for a typical aortic specimen was 5 µm, with a maximum error that was always less than 30 µm. To estimate the lowest strain error possible for the DIC configuration, the mean value of the strain was calculated for the 0 kPa case [START_REF] Reu | The art and application of DIC: Calibration sanity checks[END_REF]. The strain floor was found to be 0.6%.

Geometric Reconstruction

A deforming NURBS mesh was extracted by morphing a NURBS template to the DIC point clouds. The template was a circular domain with a diameter slightly less than that of the point cloud in the first pressure state. The NURBS surface was parameterized as a single patch containing clamped knots of 20 divisions in each parametric direction, with 22×22 control points. Since NURBS control points, in general, do not fall on the surface they describe, they cannot be directly derived from the DIC clouds. Instead, the positions of the Gauss points were obtained first using a moving least-square method [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF]. For each Gauss point, a set of nearest image points in the DIC point cloud were identified based on their distance to the Gauss point in the first pressure state. The radius of the neighboring region was automatically adjusted so that it contained at least 6 image points. The position of each Gauss point, y g , was computed using an affine interpolation

y g = y∈Ωg w j y j y∈Ωg w j (2)
where y j is the position vector for each image point in the neighborhood, Ω g , and w i is the weighting function taken to be the inverse of the distance from y j to the Gauss point. Using the same weights calculated in the first stage, the Gauss points in every pressure stage were identified.

A global least squares problem was then formulated to compute the best-fit positions of the control points. The NURBS surface was represented as

x = i N i (u 1 , u 2 )Q i (3)
where N i are the NURBS basis functions, Q i are the control points, and the pair of knot variables, (u 1 , u 2 ), represent a material point. The position of a modeled Gauss point is then given by

x g = i N i (u 1g , u 2g )Q i .
The control points were obtained by minimizing a weighted sum of x g -y g 2 over all Gauss points. This procedure was applied to each pressure state.

To assess the accuracy of the constructed NURBS mesh the distance deviation between the NURBS surface and the DIC point cloud was calculated. For each image point, y i , its closest projection point on the NURBS surface, x i , was computed [START_REF] Lu | Isogeometric contact analysis: Geometric basis and formulation for frictionless contact[END_REF]. Points whose closest projection laid outside the NURBS domain were discarded. Using the remaining points the distance deviation was calculated

d = i x i -y i 2 n 1/2 (4)
where n is the number of DIC cloud points inside the NURBS domain.

Data Analysis

Strain Surface strains were computed with the aid of the NURBS curvilinear coordinate system. The surface coordinates, u α , (α, β = 1, 2) induce a set of convected basis vectors (a 1 , a 2 ) where a α = ∂x ∂uα and x(u 1 , u 2 ) is the NURBS representation given in Eq. 3. The reciprocal basis a 1 , a 2 are computed such that a α • a β = δ β α . In the reference configuration, the basis vectors are denoted by (A 1 , A 2 ) and A 1 , A 2 .

The surface deformation gradient tensor is

F = a α ⊗ A α . (5)
It then follows that the surface Cauchy-Green deformation tensor, C, and the Green-Lagrangian strain tensor, E, are given by

C = (a α • a β ) A α ⊗ A β (6) E = 1 2 (a α • a β -A α • A β ) A α ⊗ A β . (7)
The physical components of C and E are computed by constructing a local orthonormal basis (G 1 , G 2 ) in the tangent plane spanned by (A 1 , A 2 ). The physical components of the Cauchy-Green deformation tensor, C ij , and Green-Lagrangian strain tensor, E ij , are

C ij = G i • C G j and E ij = G i • E G j , respectively.
Wall Stress For an inverse membrane boundary value problem the deformed configurations and boundary conditions are given as inputs to the finite element model and the wall stress is calculated. The balance equation that governs static equilibrium is [START_REF] Green | Large elastic deformations[END_REF][START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF])

1 √ a √ at αβ a α ,β + pn = 0 (8)
where a is det [a α • a β ], t is the Cauchy wall tension, p is the applied internal pressure, n is an outward facing unit normal, and ( ) ,β indicates ∂ ∂ u β . Note that the Cauchy wall tension t is directly related to the Cauchy stress, σ, through the current thickness of the membrane, h, via t αβ = hσ αβ = t βα . If the thickness distribution is known for every element (or indeed every Gauss point) in the mesh the local Cauchy stress can be ascertained. Using our current experimental method the thickness distribution of the tissue cannot be measured accurately. To avoid introducing errors into the stress calculation due to imprecise thickness values we have elected to use the wall tension to describe the stress state of the tissue.

The weak form of the boundary value problem reads

Ω t αβ a α • δx ,β da - Ω pn • δxda = 0 (9)
where δx is any admissible variation to the current configuration Ω. The details of the finite element procedure for solving Eq. 9 were presented in [START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF].

Briefly, the Cauchy wall tension is regarded as a function of the inverse deformation gradient. The weak form subsequently yields a set of nonlinear algebraic equations for the positions of control points in the reference configuration. While solving for the reference configuration the tension field in the current state is determined. An auxiliary material model is needed to perform the inverse analysis. The material model influences the predicted undeformed configuration; however, due to the static determinacy of Eq. 8, the influence on the wall tension is expected to be weak [START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF][START_REF] Lu | A shell-based inverse approach of stress analysis in intracranial aneurysms[END_REF][START_REF] Miller | On the prospect of patientspecific biomechanics without patient-specific properties of tissues[END_REF]. As in a previous study [START_REF] Genovese | Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder in vitro[END_REF], a neo-Hookean model was implemented. For computational efficiency, the stiffness parameter of the model was set to unrealistically high values to ensure a robust convergence.

To facilitate the inverse analysis, the control points on the outermost edge of the mesh were fixed. Since the DIC point cloud does not include the edge where the specimen was fixed to the clamp, the applied boundary conditions are only approximate. Applying any displacement based constraint in the inverse membrane analysis creates a boundary region in the solution where the stresses are inaccurate [START_REF] Lu | Inverse method of stress analysis for cerebral aneurysms[END_REF]. To minimize the influence of boundary effect, the outer ring of elements were deemed to be the boundary region and excluded from further analyses. By a retrospective comparison with the forward analysis reported later (Sec. 3.4) the size of the boundary region was confirmed.

Material Property Identification Using inverse membrane analysis, the stress was calculated at every Gauss point. Combining the stress data with the local surface strains calculated from Eqs. 5 -7, the stress-strain response at every Gauss point in the mesh is known. The local material properties at each Gauss point were then identified by fitting the local stress-strain response to a hyperelastic surface energy density. A two-dimensional version of the strain energy density proposed by Gasser, Ogden, and Holzapfel (GOH) ( 2006) was used:

w = µ 1 2 (I 1 -ln (I 2 ) -2) + µ 2 4γ e γ(I k -1) 2 -1 ( 10
)
where I 1 = tr C and I 2 = det C are the principal invariants of the Cauchy-Green deformation tensor and

I κ = C : (κI + (1 -2κ) M⊗M
) is a compound invariant consisting of isotropic and anisotropic contributions.

In the compound invariant I κ , the unit vector M = cos θ G 1 + sin θ G 2 defines the preferred orientation of the collagen fibers in the reference configuration while κ characterizes the dispersion of the collagen fibers.

The dispersion parameter, κ, varies from 0 to 1 2 . When κ = 0 all of the collagen fibers are perfectly aligned in the direction of M and when the distribution of collagen fibers has no preferential direction (isotropic), κ = 1 2 . In Eq. 10, the parameters µ 1 and µ 2 are the effective stiffnesses of the matrix and collagen fibers, respectively, both having dimensions of force per unit length. The parameter γ is a non-dimensional parameter that governs the tissue's strain stiffening response.

The second Piola-Kirchhoff wall tension, S, is written as

S = 2 ∂w ∂I 1 1 + 2 ∂w ∂I 2 I 2 C -1 + 2 ∂w ∂I κ (κI + (1 -2κ) M ⊗ M) . ( 11
)
Substituting Eq. 10 into Eq. 11 one finds

S = µ 1 I -C -1 + µ 2 e γ(Iκ-1) 2 (I κ -1) (κI + (1 -2κ) M ⊗ M) (12)
noting that the second Piola-Kirchoff wall tension is related to the Cauchy wall tension via t = 1 √ I2 F S F T . The values of the model parameters µ 1 , µ 2 , γ, κ, and θ were determined by minimizing the sum of the squares difference between the stress computed from the inverse membrane analysis and those computed using Eq. 12. The nonlinear minimization was solved in MATLAB (MathWorks,v. 7.14) where the model parameters were constrained such that: µ 1 , µ 2 , γ > 0, 0 ≤ θ ≤ π, and 0 ≤ κ ≤ 1 2 . To ensure that the minimization routine converged to a global minimum, three sets of initial guesses were used. Since the parameters θ and κ are tightly bound, the upper bound, lower bound, and middle value were used. For the parameters µ 1 , µ 2 , and γ random values between 0.01 and 5 were generated. The parameter identification was relatively unaffected by changes to the initial guess and consistently converged to the same solution.

Forward Validation

The pointwise wall tensions and property distributions were validated using forward finite element analysis. This was accomplished by mapping the pointwise mechanical properties back to the NURBS mesh and simulating the inflation process. The configuration at a pressure of 3 kPa was taken as the reference configuration. The applied pressure was assumed to balance the weight of the tissue in this configuration rendering it stress-free. Since the material parameters in the perimeter ring of elements were not identified, these elements were excluded from the forward analysis. Displacement boundary conditions, derived from the experimental DIC data, were applied to these elements. The wall tensions, strains, and displacements calculated in the forward analysis were compared to the DIC displacements, the DIC strains computed on the NURBS mesh, and the stresses from the inverse analysis. The relative error was computed at every Gauss point and averaged over the computational domain. To illustrate the influence of assuming pointwise heterogeneous properties, two additional forward validations were performed using an elementwise and a homogeneous distribution of mechanical properties. To identify the elementwise material properties the nonlinear minimization described earlier was applied to each of the interior elements. For the homogeneous case, all of the pointwise material properties were averaged to identify one set of material properties. The forward validations for the elementwise and the homogeneous cases were performed using the same method described for the pointwise case.

Results

A bulge inflation test to failure was performed on an ATAA collected from a male patient who was 55 years old. The diameter of the aneurysm as determined by pre-surgical CT scan was 55 mm. The sample ruptured at 117 kPa. However, it must be noted that, in terms of the stress it induces the pressure applied in the test is not equivalent to a blood pressure, for example, 15 kPa measured by the pressure sensor is not equivalent to a blood pressure of 15 kPa (113 mmHg) due to the new shape of the specimen. The pressure and DIC data collected during the bulge inflation test were used to generate a deforming NURBS mesh and identify the local stress-strain response during the inflation test. Using the pointwise stress-strain data, the spatial distribution of the mechanical properties was identified.

Geometric Reconstruction

Using the experimental DIC point cloud a deforming NURBS mesh was generated of the ATAA sample. Figure 2 shows the extracted NURBS mesh at a pressure of 117 kPa superimposed on the experimental DIC point cloud at the same pressure. Using Eq. 4 the distance deviation for each pressure stage was computed. As shown in Fig. 2, the overall the distance deviations are very small, however the values are significantly larger at low pressures. The distance deviation rapidly decreases as the tissue is inflated and stabilizes by 30 kPa.

Local Stress and Strain Response

Figure 3 shows the distributions of the magnitude of the Cauchy wall tension, t, and Green-Lagrangian strain, E, at an applied pressure of 30 kPa (Fig. 3a,b) and 75 kPa (Fig. 3c,d). The distribution of wall tension and strain remained similar throughout the inflation of the specimen. Components of the Cauchy wall tension and Green Lagrange strain at the center of each element in the computational domain are provided in an on-line supplementary data file.

In general, at each Gauss point both the normal strains and the planar shear strains were non-zero. To facilitate plotting of the local stress-strain response, the axes of principal strain were identified and the local stresses and strains were rotated into the principal strain axes. In Fig. 4, the three components of the wall tension in the principal strain axes, t11 , t12 , and t22 are plotted against the principal stretches λ 1 and λ 2 . As expected, the local stress-strain response shows the non-linear stiffening behavior that is common in arteries. In a small region where rupture eventually occurred the ATAA appears to yield (Fig. 4c). The locations of this localized yielding correspond to strain concentrations in zone where rupture initiates (Fig. 7c).

Material Property Identification

The proposed model for the elastic behavior of the ATAA was able to fit the bulge inflation data well (0.81 < R 2 < 0.99). Lower values of the correlation coefficient were located in the small zone where rupture eventually occurred. Excluding this region the minimum value of R 2 was 0.96. The experimental data (points) and model fits (lines) for three Gauss points are shown in Fig. 4. The distributions of the material parameters are plotted in Fig. 5 and they clearly display a heterogeneous distribution. The parameter µ 1 displayed the sharpest changes in value while the parameters µ 2 , κ, and γ changed more gradually. Not surprisingly, the values of µ 2 are an order of magnitude larger than µ 1 reflecting the difference in stiffness between the collagen fibers and matrix. The values of κ are approximately 0.5 in the center indicating that the tissue as a whole is only moderately anisotropic. The non-uniform distribution of κ suggests that there are regional differences in the fiber organization. In Fig. 5e, the angle θ that defines the primary orientation of the collagen fibers is plotted. Note that θ is defined locally relative to the local basis vector G 1 which is parallel to the horizontal meshlines. Keep in mind that when the value of κ is approximately 0.5, the fiber orientation θ is irrelevant since there is no preferred fiber direction. The algorithmic mean and standard deviation for each of the material parameters are recorded in Table 1.

Forward Validation

Using the pointwise material properties excellent agreement was achieved between the forward finite element prediction and the experimental results. In Fig. 6, the predicted and reconstructed geometries at four pressures are presented. While there are visible deviations at low pressure (Fig. 6a), the configurations at higher pressures are indistinguishable (Fig. 6b-d). The predicted magnitude of the wall tension and strain were also very similar to their corresponding experimental results, as illustrated in Fig. 7. The quality of the forward prediction is summarized in Table 2 where the average percent error for the wall tension, strain, and displacement are presented at 8 pressures. The highest errors between the forward predictions and experimental data occur at low pressures. The errors then rapidly decline, stabilizing at approximately 60 kPa reflecting the good agreement seen in Figs. 6 and7. As expected, the smallest errors were found in the displacement predictions. These errors in the displacement then lead to the larger errors observed in the wall tension and strain. Similar results were found using elementwise material properties as shown in Table 3. Although the errors in displacement, wall tension, and strain are elevated compared to the pointwise case, the mean error for displacement, wall tension, and strain were all below 6% when the error stabilized. The small magnitude of the errors confirm that the material property distributions were accurately identified.

When the homogeneous material properties given in Table 1 were used, the forward simulation produced large displacement errors, which in turn, generated large errors in the strain and wall tension. The average displacement error for the homogeneous case was approximately 10% when the error stabilized and the resulting average strain and wall tension errors were on the order of 20%. The large errors that result from assuming homogeneous material properties suggest that the spatial variation in the material parameters must be accounted for in order to accurately replicate the behavior of the ATAA.

Discussion

In this manuscript, we have presented a novel method to identify the pointwise heterogeneous distribution of material properties in soft tissues. The feasibility of the method was demonstrated using a human ATAA. However, in practice, this method could be used to characterize a variety of thin soft tissues such as skin, aortic valves, or the bladder. By combining DIC, bulge inflation testing, and inverse membrane analysis we are able to obtain the stress-strain response locally, at every Gauss point. Using this local stress-strain data the elastic properties were identified. Specifically for the ATAA a modified form of the GOH model was fitted to the local stress-strain curves to identify the heterogeneous distribution of the material properties. A forward validation was also performed to show that a heterogeneous distribution of the material properties must be used to accurately describe the response of the ATAA.

In previous studies using bulge inflation tests, the membrane stresses were frequently identified using the Laplace equation which requires that the membrane is axisymmetric. Investigators assumed that the tissue deformed into either a spherical cap [START_REF] Mohan | Failure properties of passive human aortic tissue. II-Biaxial tension tests[END_REF] or an elliptical cap [START_REF] Coudrillier | Biomechanics of the human posterior sclera: Age and glaucoma related changes measured using inflation testing[END_REF][START_REF] Marra | Elastic and rupture properties of porcine aortic tissue measured using inflation testing[END_REF] so that the Laplace equation could be used. However, during our bulge inflation tests, the ATAA did not maintain an axisymmetric shape but, instead, deformed into a prolate dome without any pronounced symmetry (Fig. 6d). For membranes with no marked symmetry, the Laplace equations cannot provide accurate estimations of the stress and the numerical approach presented in this manuscript is necessary.

The protocol presented, in particular the inverse membrane analysis, is quite sensitive to the quality of the displacement data [START_REF] Zhao | Pointwise identification of elastic properties in nonlinear hyperelastic membranes -Part II: Experimental validation[END_REF]. The accuracy of the stress reconstruction depends heavily on local changes in the surface curvature which are calculated from the DIC displacement data. Numerical derivatives are performed on the displacement data to calculate the strain so, any noise in the displacement data is transferred to and amplified in the strain data. When the surface is non-smooth, having small local undulations, the accuracy of the stress and strain reconstruction is also compromised. To overcome these possible accuracy issues, we have employed a stereoscopic DIC system where we have estimated that the displacement measurements are accurate to 30 µm (Eq. 1) and second order NURBS elements. The NURBS approach uses isogeometric analysis and results in an intrinsically smooth method to describe the geometry [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], which, in turn, produces smooth stress and strain data. In Fig. 2, the higher distance deviations in the geometric reconstruction (Eq. 4) at pressures below 30 kPa are believed occur due to a non-smooth surface geometry. As the ATAA is inflated, the specimen surface becomes smooth. Reflecting this change in the geometry the mean distance deviation gradually declines and stabilizes around 30 kPa. The results of the forward validation (Figs. 67) confirm that the displacement, strain, and stress were correctly identified using our protocol.

In general the quality of the fits performed using the GOH model were very good but at the low pressures the fits were sometimes poor as seen in Fig. 4. Several factors contributed to the decreased accuracy at low pressures. First, the assumption that the configuration at 3 kPa is stress-free is debatable. The flattening of the ATAA sample from its cylindrical state and the inversion so that the lumen faced out could conceivably lead the assumed reference state to be stressed. Second, by defining the objective function as the sum of the squares difference between the stresses computed from the inverse simulation and the model stresses automatically assigns equal weight to each state, in effect skewing the fits toward the high stress terms. Consequently poor fits at low values of stress and strain do not significantly change the overall value of the objective function. The introduction of a non-uniformly weighted objective function could improve the fits by assigning higher weights to stress-strain points at low pressures.

The material properties identified for the ATAA are highly heterogeneous as evidenced by Fig. 5. The values of the parameters µ 1 , µ 2 , γ, and κ ranged from 0.002-0.04 N/mm, 0.15-0.65 N/mm, 1.00-2.50, and 0.30-0.46, respectively. The map of the dispersion parameter, κ, in Fig. 5d suggests that the ATAA is regionally anisotropic. This finding may explain the conflicting reports in the literature on the whether the thoracic aorta shows directional dependence [START_REF] Choudhury | Local mechanical and structural properties of healthy and diseased human ascending aorta tissue[END_REF][START_REF] Okamoto | Mechanical properties of dilated human ascending aorta[END_REF][START_REF] Pham | Biomechanical characterization of ascending aortic aneurysm with concomitant bicuspid aortic valve and bovine aortic arch[END_REF]Zhou and Fung, 1997). Figure 5d reveals that the degree of anisotropy observed is location dependent. This observation could be confirmed by studying of the microstructural organization of the ATAA using histology to independently quantify the values of the structural parameters κ and θ.

To show that it is necessary to consider the heterogeneous distribution of the mechanical properties, we performed three different forward finite element simulations using pointwise, elementwise, and homogeneous material properties. In Tables 23, the errors in displacement, wall tension, and strain for the pointwise and elementwise cases are given. It is clear that homogenization leads to increased errors in the forward finite element simulations. While the errors in the final stage for the elementwise case were 5 times greater than for the pointwise case, the mean errors for displacement, strain, and wall tension were only 0.84±0.93%, 4.72±5.72%, and 5.49±5.70% respectively. Note that in this case, the use of elementwise properties still requires that the heterogeneities of the material parameters are known at a spatial resolution of 1 mm. For the homogeneous case the mean displacement error by the final stage had only decreased to 7.49%. Large displacement errors (> 5%) compromise the stress calculations as discussed earlier. Our results suggest that assumption of homogeneity for the ATAA may inhibit researchers' ability to accurately recreate physiologic stresses and strains in finite element simulations.

Although the primary aim of this method is to identify the local elastic response of soft tissues, it is also an apt method for the identification of rupture properties. For our study case of the ATAA, the most commonly used method to estimate the rupture stress is uniaxial tensile testing [START_REF] Iliopoulos | Regional and directional variations in the mechanical properties of ascending thoracic aortic aneurysms[END_REF][START_REF] Sokolis | Effect of layer heterogeneity on the biomechanical properties of ascending thoracic aortic aneurysms[END_REF][START_REF] Vorp | Effect of aneurysm on the tensile strength and biomechanical behavior of the ascending thoracic aorta[END_REF]). Yet it is well known that the in vivo loading experienced by arterial tissues is much better approximated by biaxial tests [START_REF] Marra | Elastic and rupture properties of porcine aortic tissue measured using inflation testing[END_REF]. By design, our approach can easily identify the rupture stress in a non-uniform two-dimensional stress environment, and therefore, the rupture stresses calculated are believed to be more accurate than those identified from uniaxial tests. In addition, our method leaves open the possibility of evaluating the mode of rupture by characterizing how the initial tear propagates in the tissue.

There are, however, several limitations to our method. First, this protocol is only valid for membranes. From the standpoint of mechanics, a membrane is a solid structure where the thickness is small compared to the other dimensions and the structure's resistance to bending is small. At low pressures (p < 20 kPa), the membrane condition may be compromised due to a relatively small length to thickness ratio when the specimen is only slightly inflated. It is also possible that the bending effects in the ATAA may not be negligible at low pressures where bending is dominant mode of deformation. The nonlinear elastic behavior of the ATAA may amplify the role bending plays in the stress-strain response. Indeed, this may be one of the causes of the large errors that were found in the forward validation at low pressures. Second, only a limited variety of deformations can be produced using a bulge inflation test. For a bulge inflation test one does not prescribe a biaxial deformation or stress state instead, a pressure is specified and the resulting deformation recorded. Therefore it is unclear if the parameters from the bulge inflation test can be used to predict the response of the ATAA in uniaxial tension, for example. Finally, we only considered three types of material property distributions: pointwise, elementwise, and homogeneous. To truly identify the spatial resolution necessary to accurately describe the heterogeneous properties of the ATAA other levels of refinement should be considered. Despite these limitations, our approach provides a novel method that can characterize the complex heterogeneity observed in soft tissues.

In conclusion, we have developed a protocol for characterizing the heterogeneous mechanical properties of thin planar tissues. The true strength of this combined experimental and computational approach is that the method is unaffected by the complexity of the heterogeneity or the nonlinear elastic behavior of the sample. Although the method was applied to the ATAA in the feasibility study, the method could be widely applied to a variety of planar soft tissues. In the future, we plan to conduct a large population study on the ATAA. This study would combine our method with local structural analysis such as histology or multi-photon microscopy to relate the microstructural organization of the tissue to the distribution of the mechanical properties. The additional structural data is necessary to move from merely measuring the heterogeneous nature of the material properties to uncovering the links that exist between the measured properties and the underlying microstructure. We also plan to use the heterogeneous properties from the population study to quantify the influence that material heterogeneity has on the predicted peak wall stress in patient specific finite element models. Furthermore, this method could be transferred to a clinical setting in the future to estimate the material properties of soft tissues in vivo. When coupled with an appropriate imaging modality (e.g. 4D MRI) the evolution of the material properties as the ATAA enlarges could be measured on a patient-to-patient basis as a part of the surveillance routine. 10(2):177-189. Zhou J, Fung YC (1997) The degree of nonlinearity and anisotropy of blood vessel elasticity.

Proc Natl Acad Sci, 94(26):14255-14260. doi:10.1073/pnas.94.26.14255

Fig. 1

 1 Fig. 1 Experimental setup and test sample (a) before testing and (b) after rupture

Fig. 2

 2 Fig.2NURBS mesh (blue) superimposed on the experimental DIC point cloud at a pressure of 117 kPa along with a plot of the distance deviation calculated using Eq. 4 as a function of pressure.

Fig. 3

 3 Fig. 3 Contours of the magnitude of the wall tension (N/mm) and Green-Lagrange strain at (a,b) 30 kPa and (c,d) 75 kPa.

Fig. 4

 4 Fig.4Local wall tension versus stretch curves at three different locations. Note that the curves have been rotated into their principal strain state to facilitate plotting.

Fig. 5

 5 Fig. 5 Distribution of the identified material parameters over the ATAA.

Fig. 6

 6 Fig.6Comparison of the geometry constructed from the DIC point clouds (blue) and those predicted from the forward finite element analysis using the pointwise material properties (yellow) at (a) 15 kPa, (b) 30 kPa, (c) 75 kPa, and (d) 117 kPa. Note that the geometry is almost a perfect match leading the blue and yellow lines to overlap.

Fig. 7

 7 Fig.7Comparison of the magnitude of the (a,b) wall tension (N/mm) and (c,d) Green-Lagrange strain from (a,c) the experimental data and (b,d) the forward finite element prediction at a pressure of 117 kPa. Note that the red ellipse in the top center of (c) is the region where rupture initiated.

Table 1

 1 Mean value and standard deviation of the material parameters.

	Material Parameters	Mean ± Std
	µ 1 (N/mm)	0.0122 ± 0.0111
	µ 2 (N/mm)	0.2747 ± 0.1051
	γ	1.6090 ± 0.3527
	κ	0.4019 ± 0.0478
	θ (rad)	1.5560 ± 0.6003

Table 2

 2 Comparison of forward finite element predictions using pointwise material properties with the experimental measurements

	Pressure (kPa)	Displacement	Wall Tension		Strain
		Max. (mm) Error (%)	Max. (N/mm) Error (%)	Max. Error (%)
	15	7.15	1.94±1.68	0.21	8.25±6.54	0.39	7.99±6.33
	30	9.63	0.42±0.42	0.39	3.92±3.35	0.61	2.67±2.50
	45	10.75	0.27±0.25	0.61	3.30±2.84	0.71	1.50±1.44
	60	11.49	0.13±0.14	0.84	2.78±2.45	0.77	0.95±1.00
	70	12.10	0.10±0.11	1.08	2.18±2.00	0.81	0.73±0.85
	90	12.62	0.09±0.11	1.33	1.69±1.68	0.85	0.69±0.85
	105	13.15	0.12±0.14	1.59	1.61±1.73	0.88	0.79±1.08
	117	14.11	0.11±0.17	1.86	2.22±1.95	0.94	0.91±1.45

Table 3

 3 Comparison of forward finite element predictions using elementwise material properties with the experimental measurements

	Pressure (kPa)	Displacement	Wall Tension		Strain
		Max. (mm) Error (%)	Max. (N/mm)	Error (%)	Max.	Error (%)
	15	7.14	2.41±2.68	0.21	9.69±10.27	0.39	9.30±10.06
	30	9.63	1.89±2.04	0.39	7.52±7.44	0.61	6.30±6.74
	45	10.75	1.03±1.09	0.61	6.14±5.96	0.71	4.58±5.29
	60	11.49	0.51±0.59	0.84	5.43±5.62	0.77	3.85±4.93
	70	12.10	0.48±0.59	1.08	5.11±5.83	0.81	3.86±4.97
	90	12.62	0.67±0.77	1.33	5.02±6.02	0.85	4.17±5.19
	105	13.15	0.85±0.95	1.59	5.56±6.34	0.88	4.61±5.55
	117	14.11	0.85±0.93	1.86	5.49±5.70	0.94	4.72±5.72

Acknowledgements This work was supported in part by The University of Iowa Mathematical and Physical Sciences Funding Program. Dr. Davis was supported by the Whitaker International Scholars Program. The authors thank Dr. Aaron Romo for his assistance with the experimental protocol.


José Manuel Sánchez

Imen Grida

Ben Yahia

Noel Crespi

Self-Modeling Based Diagnosis of Software-Defined Networks

Keywords: self-modeling, self-diagnosis, Bayesian networks, SDN, NFV, Case1: 1.9 Case2: 1.8 Case1: 1 Case2: 1.1 Case1: 1 Case2: 1.1 Case1: 0 Case2: 31.1 Case1: 1.9 Case2: 1.8

Networks built using SDN (Software-Defined Networks) and NFV (Network Functions Virtualization) approaches are expected to face several challenges such as scalability, robustness and resiliency. In this paper, we propose a self-modeling based diagnosis to enable resilient networks in the context of SDN and NFV. We focus on solving two major problems: On the one hand, we lack today of a model or template that describes the managed elements in the context of SDN and NFV. On the other hand, the highly dynamic networks enabled by the softwarisation require the generation at runtime of a diagnosis model from which the root causes can be identified. In this paper, we propose finer granular templates that do not only model network nodes but also their sub-components for a more detailed diagnosis suitable in the SDN and NFV context. In addition, we specify and validate a self-modeling based diagnosis using Bayesian Networks. This approach differs from the state of the art in the discovery of network and service dependencies at run-time and the building of the diagnosis model of any SDN infrastructure using our templates.

I. INTRODUCTION

SDN (Software-Defined Networks) and NFV (Network Function Virtualization) is a novel phenomenon that is on the wish list of major industrial players (vendors, operators, content providers, software editors) as the means to achieve greater flexibility in managing the network, faster service deployment and provisioning while reducing operational costs. SDN is expected to pave the way towards network programmability by proposing network architecture based on abstraction, open interfaces, and control plane-data plane separation. On the other hand, NFV promise is to turn traditional network functions into virtual ones called Virtual Network Functions (VNF) and embed them into commoditized hardware. The combination of both approaches is commonly agreed to be the best solution and is the next industry move despite the preliminary stage of the management of such environment [START_REF] Esteves | On the management of virtual networks[END_REF] [START_REF] Kreutz | Software-Defined Networking: A Comprehensive Survey[END_REF]. In particular, resiliency properties become fundamental for both technologies as discussed in [START_REF] Sanchez | Softwarized 5G networks resiliency with self-healing[END_REF]. We rely on self-diagnosis [START_REF] Kephart | The Vision of Autonomic Computing[END_REF] as a first step towards resiliency. Self-diagnosis is an autonomic capability where the network is aware of any abnormal state and diagnoses itself to determine the root cause to perform the appropriate remediation or recovery actions. Those actions may be based on redundancy or diversity mechanisms, which bypass the presumed faulty network elements.

In the SDN and NFV context, the highly changing networks impose numerous diagnosis challenges, specifically how to detect continuous changes and update dependencies among virtual and physical resources. The contribution of this paper is twofold: 1) definition of finer granularity templates that model the dependencies among SDN nodes (physical and logical) as well as smaller subcomponents inside them (e.g. CPU, network cards, etc).

2) Specification and validation of a self-modeling approach that tracks changes on the network topology of the networked nodes and their corresponding VNFs at runtime. The structure of the paper is as follows: Section II motivates our work and Section III summarizes the related work on model-based network diagnosis which includes the diagnosis techniques and the models that are in use. Section IV presents the proposed self-model based diagnosis framework including the template definition and the defined algorithm. Section V presents the experimental validation. Finally, Section VI summarizes our findings and outlines the future work.

II. MOTIVATION

We consider an end-to-end service delivered to clients through SDN and NFV based networks. This combination of SDN and NFV raises the following questions with respect to diagnosis:  The first challenge is how to model the dependencies among VNFs at runtime, as these dependencies may vary over time and depend on the service contracted by each user (Fig. 1 (a)).  The second challenge is how to model at runtime the dependencies of the underlying network topology, when it is considered dynamic (Fig. 1 (b)). The network topology is dynamic due to several reasons like the connections and disconnections of users to the Access Points (AP), handovers, and especially VNF migrations and the orchestration of new services.

978-1-4799-7899-1/15/$31.00 ©2015 IEEE  The third challenge is how to model the two types of control in SDN, out-of-band and in-band. In out-of-band control, the controller is directly connected to every switch by a dedicated control link, but in in-band control, the controller is only directly connected to the master switch, becoming this one an intermediate node that connects the controller with the rest of switches (slaves). As it will be explained in section III, the most popular approach for diagnosis is the model-based diagnosis. Our approach targets to answer the challenges of diagnosis in the context of SDN and NFV and tracks the dynamic network topology. For this, we propose 1) a model that is suitable to handle the specific elements to monitor within a combined SDN and NFV scenario, and 2) a self-modeling algorithm to enable the automation of the diagnosis model.

III. RELATED WORK ON MODEL-BASED DIAGNOSIS

The diagnosis model is generally a dependency graph. In general, the dependency graph is manually generated from an operational team"s knowledge. This manual generation is valid for static network topologies, but not for dynamic and elastic networks such as those expected with SDN and NFV. A selfmodeling approach is presented in the literature as the automatic generation of this model [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF]. The dependency graph is first generated from a given data set (databases [START_REF] Singh | Automated probabilistic modelling for relational data[END_REF], genetic algorithms [START_REF] Kabli | A Chain-Model Genetic Algorithm for Bayesian Network Structure Learning[END_REF], or ontologies [START_REF] Fenz | An ontology-and bayesian-based approach for determining threat probabilities[END_REF]). An inference engine (based on algorithms such as Bayesian Networks [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF][10], Occam"s Razor [START_REF] Bahl | Towards highly reliable enterprise network services via inference of multi-level dependencies[END_REF], or others) reasons then the dependency graph in order to extract the root cause.

In this paper, we focus on the generation at runtime of the dependency graph from the dynamic network topology. We consider the following works [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF][10] [START_REF] Kandula | Detailed diagnosis in enterprise networks[END_REF][12] which we explain and compare hereafter:

Hounkonnou et al. in [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF] propose a self-modeling approach based on patterns to enhance Bayesian Networks (BN) and apply it to diagnose the IP Multimedia Subsystem (IMS). This approach generates offline a generic model (pattern) that is based on the four IMS layers (physical, functional, procedural and service). A pattern describes the dependencies among resources used by an IMS service. When a failure occurs in the IMS service, the algorithm locates online instances of that pattern in a given network topology and it generates the corresponding BN instance. However, this work assumes that the network topology remains static during the diagnosis process. Furthermore, the granularity of the diagnosis covers the four layers of IMS but it diagnoses only the network resource level, without considering smaller sub-components inside them. It is also worth mentioning that it does not consider virtual environments.

Bennacer et al. in [START_REF] Bennacer | Optimization of fault diagnosis based on the combination of Bayesian Networks and Case-Based Reasoning[END_REF] base their self-modeling approach on Chi-squared statistical tests. These tests learn the dependencies among modeled variables. They utilize a "significance level" to decide their dependencies. This selfmodeling approach considers physical symptoms as variables on each network node. However, the granularity of the diagnosis remains at the network node level, where smaller sub-components are not considered. Furthermore, the diagnosis is focused on the physical network nodes and is not considering the logical resources (e.g. virtual resources) running over them. In addition, inappropriate values of "significance level" may lead to errors when building the dependency graph. Kandula et. al. in [11] present a self-modeling approach based on templates. They define specific templates for each element such as a machine, an application process, a neighbor set and a path. Each template is characterized by several state variables to achieve a detailed diagnosis. Same as [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF][10], the model granularity remains at the node level and does not consider smaller sub-components inside them. Bahl et. al. in [START_REF] Bahl | Towards highly reliable enterprise network services via inference of multi-level dependencies[END_REF] models the dependencies among different services, but also the dependencies among network nodes given by the network topology. However, they do not consider smaller sub-components inside nodes or virtual resources. This article advances the state of the art by describing a selfmodeling based diagnosis to discover at runtime the dependency model of SDN-NFV infrastructures. It then models automatically an SDN based end-to-end service, its underlying network topology, and the type of control implemented (in-band or out-of-band). Compared to [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF], which assumes that the network topology, services and configuration remain static, our proposal, in the context of SDN-NFV, assumes a continuous changing network topology. We then create the dependency graph by instantiating templates at runtime. Unlike [START_REF] Bennacer | Optimization of fault diagnosis based on the combination of Bayesian Networks and Case-Based Reasoning[END_REF], our self-modeling approach builds the dependency graph from the network topology, instead of statistical tests to avoid inaccuracies in the model. The idea of using different templates to describe different network elements was inspired by [START_REF] Kandula | Detailed diagnosis in enterprise networks[END_REF]. Our templates are tailored to suit the SDN elements" particularities and VNFs deployed in the network nodes. Contrary to [START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF][10] [START_REF] Kandula | Detailed diagnosis in enterprise networks[END_REF][12], we consider a finer dependency model granularity, that considers smaller sub-components within a network node and their internal dependencies (physical and logical) to be able to diagnose a combined SDN and NFV environment where virtual and physical parts such as network cards and CPU need to be diagnosed. We use Bayesian Networks approach [START_REF] Ben-Gal | Bayesian Networks[END_REF] for the root cause calculation which is enhanced with our selfmodeling approach.

IV. SELF-MODELING BASED DIAGNOSIS FRAMEWORK

We present a self-modeling based diagnosis framework for SDN (Fig. 2). We propose: (1) a template to model SDN end-to end service

A P 2 A P 1 UE 1 UE 1 UAP 1 UAP 1 VNF 1 VNF 1 VNF 3 VNF 3 VNF 2 VNF 2 UAP 1 UA 1 VNF 3 VNF 1 H 1 H 2 H 3 H 4 H 5 H 6 H 7 R 1 R 2 R 3 R 4 R 6 R 7 controller User: U 1 VNF 2 UA 1 control plane data plane (a) (b)
elements and (2) a self-modeling approach that builds automatically the model from the network topology. The self-modeling building block is composed of the topology interpreter (a) and the dependency graph building (b) subblocks. This block automatically builds the dependency graph based on the network topology and the type of SDN control (in-band and out-of-band).

The root cause calculation building block finds the root cause through BNs. It receives a service alarm about service degradations or unavailability and correlates this alarm with network observations (the state of the physical and logical resources of each network node).

A. Background on Bayesian Networks

BNs utilize probabilistic properties to perform reasoning in uncertain domains [START_REF] Ben-Gal | Bayesian Networks[END_REF]. The model, the dependency graph, is represented by a set of vertices V describing events interlinked by edges that represent the dependencies among vertices. The pair describes a BN. is the dependency graph, and contains the parameters of the BN, which take shape of Conditional Probability Tables (CPT) that specify the probability of every child vertex of the dependency graph given all value combinations of its parents. The prior probability of failure (p) is different for each vertex. The dependency graph G is composed of observable and nonobservable vertices. To reason over the dependency graph, we set the network observations into its observable vertices as evidence and the BN algorithm determines the most likely root cause(s).

B. Templates for modeling SDN networked elements

We define a network element as any type of network nodes and links. We propose a template for each network element, so a template for network nodes and another for links. These templates describe the relationships between virtual and physical sub-components inside each network element.

1) Network node template The template of a network node (TN n in Fig. 3(a)) is composed of a physical layer and a logical layer, following the TMF Information framework specifications [START_REF]TMF Information Framework[END_REF].  The physical layer encompasses the state of physical resources such as CPU and network cards. We consider two states for those physical sub-components inside network nodes: up or down.  The logical layer encompasses the state of the VNFs running inside each node. We consider that network nodes perform one or several VNFs. For instance, the controller runs the appropriate network function that installs the rules in the switches as a VNF. A VNF relies on a software process, with suitable configuration settings. The logical layer contains three sub-layers in accordance with a threestate life-cycle for each VNF: initiated, configured, and activated. Initiated: where the underlying software process of the network function is launched; configured: where the VNF is set with the optimal attributes to perform the network function; and activated: where the network function is ready to accept any request. The layers of the templates are manually predefined, but the number of network cards or running VNF per network node are retrieved from the network topology.

(controllers (C 1 … C n ), slave switches (SS 1 … SS n ), master switches (MS 1 … MS n ), hosts (H 1 … H n )) is an instance of this template.
2) Link template The template of a network link (TL n Fig. 3 (b)) is simpler and it is composed of the physical layer and one vertex. Each type of link (control links (CL 1 … CL n ), access links (AL 1 … AL n ), and inter switch links (IL 1 … IL n )) is an instance of this template.

C. Self-modeling building block

The self-modeling building block automatically models outof-band and in-band SDN networks of any topology, by parsing the input data that contains the network topology, and automating the model creation process.

1) Topology Interpreter: The topology interpreter is a northbound application that requests the network topology from the controller. The controller answers this request by providing the network topology in a JSON format. The topology data structure depends on the controller: for example, OpenDaylight and a Floodlight controller provide two different data structures, differing in the number and type of fields and the field names (Fig. 4). {"hostConfig": [{"dataLayerAddress":"00:00:00:00:00:01", "nodeType":"OF", "nodeId":"00:00:00:00:00:00:00:01", "nodeConnectorType":"OF", "nodeConnectorId":"1", "vlan":"0","staticHost":false, "networkAddress":"10.0.0.1"}]}

[{"entityClass":"DefaultEntityClass", "mac":["00:00:00:00:00:01"], "ipv4":["10.0.0.1"],"vlan":[], "attachmentPoint": [{"switchDPID":"00:00:00:00:00:00:00:01" , "port":1,"errorStatus":null}], "lastSeen":1401877225763}]

Fig. 4. JSON data structures provided by: (a) OpenDaylight, (b) Floodlight Thus, in a second stage, the topology interpreter classifies each network element into one of the following types (controllers (C j ), slave switches (SS j ), master switches (MS j ), hosts (H j ), control links (CL j ), access links (AL j ), and inter switch links (IL j ). The result is then the network descriptor, which contains the network elements (nodes and links) at instant t, and the link descriptor, which specifies the end points of each link.

2) Dependency graph building block: This block builds the global dependency graph (G) from the network descriptor. It is based on a three-step algorithm which we call: Template instantiation, topological sorting, and edge addition.

I) The template instantiation algorithm This algorithm receives as input the network descriptor. It provides as output the global dependency graph G composed of the dependency graphs (GL n or GN n ) built from the templates of each network element. It follows this methodology for each network element found in the network descriptor:  Identifies the type of network element  Instantiates its corresponding template (TLn for links and TNn for nodes)  Builds its associated dependency graph (GLn for links and GNn for nodes) from the instantiated template  Appends the dependency graph into the global dependency graph G. Each dependency graph (GL n or GN n ) contains edges called here E INSIDE . E INSIDE edges (in dash in Fig. 3) depict the dependencies among sub-components inside each dependency graph (GL n or GN n ). Each vertex in GL n or GN n is a state variable described by a different CPT. A probability of failure (P FAILURE ) is defined for each vertex even though the parents vertices are functioning. As an example, P FAILURE for the network card (NC) vertex, P(NC=down) is p, despite the proper functioning of its parent vertex (vCPU). The vertices of the global dependency graph G are not topologically sorted yet. We call these vertices as V UNSORTED .

II) The topological sorting algorithm

It sorts topologically the vertices of the global dependency graph G. It receives as input the dependency graph G with non-topologically ordered vertices (V UNSORTED ) and it provides as output the same global dependency graph but topologically sorted (V SORTED ). As an example, GN n and GL n (Fig. 3) present a topological order when separated, because any vertex index is repeated and parents are numbered before children vertices. However, when both are combined to obtain the global dependency graph G, the topological order is not respected, as both GN n and GL n contain repeated vertex indexes (e.g. both contain value "1"). The topological sorting algorithm solves this issue.

V. VALIDATION OF DIAGNOSIS MODULE

We test our self-modeling based diagnosis in a centralized SDN architecture based on a Floodlight controller. This module runs over the controller for two reasons: (1) to have a global view of the network, and (2) to keep the diagnosis framework independent from any specific southbound protocol. The network topology is obtained via the northbound interface (REST API) trough passive monitoring, to avoid introducing traffic overhead like ping tool. We use Mininet to simulate the SDN network. First, we prove that our self-modeling algorithm can interpret both the topology and the control type of SDN (out-of-band and in-band). Next, we study the scalability of this algorithm and finally we validate the diagnosis results and their variation under changing network conditions.

A. Self-modeling Validation

We test the model generation of a linear topology with two switches and a controller with two hosts connected under outof-band (Fig. 5 (a)) and in-band control (Fig. 5 (b)). in-band Fig. 5 shows the resulting global dependency graph G built by the self-modeling algorithm. This graph G is topologically ordered. We analyze the output of the algorithm hereafter:

-The self-modeling algorithm interprets the type of control:

In out-of-band control, the self-modeling algorithm instantiates two control links (instances: CL 1 , CL 2 , vertices: 2, 3), where CL 1 connects both the network card of the master switch (instance: MS 1 , vertex: 13) and the network card of the controller (instance: C 1 , vertex: 11). The controller template (C 1 ) has two network cards connected to two switches (vertices: 11, 12) through the control link instances CL 1 and CL 2 . In in-band control, it only instances one control link (instance: CL 1 , vertex: 2) because the controller is only connected to the master switch (MS 1 ). The controller template then has one network card (vertex: 10), which is connected to the network card of the master switch instance (vertex: 11) through the control link instance. The other switch is slave (instance: SS 1 ) and communicates to the controller trough the link IL 1 that connects to the master switch.

-The self-modeling algorithm interprets the network topology: For both types for control, it connects both switches MS 1 and MS 2 through the inter switch link instance (IL 1 We study the speed of the self-modeling algorithm as a function of the number of network elements (N ELEMENTS ). We launched the self-modeling algorithm for both linear and tree topologies in out-of-band control, ranging from 15 up to 500 network elements. We averaged the computing time 20 times per topology to obtain more reliable figures. Fig. 6 shows an exponential trend in the growth of self-modeling time with the number of elements for both topologies. Linear topologies scale better than tree topologies, but in both cases the selfmodeling time remains less than 30 seconds.

B. Validation of root cause analysis

We validate the diagnosis module in a linear topology (Fig. 7) in out-of-band control in two different scenarios.

In the first scenario, we study two cases, where we force certain failures on the SDN infrastructure and test if the BN engine can pinpoint accurately the failed element. It receives an alarm about this failure on SDN infrastructure and it is fed with the states of network cards of all the network nodes as network observations. The BN engine correlates these observations with this alarm to find the root cause.

In case 1 (Fig. 7(a)), the actual root cause is a total shutdown of the controller. The BN engine determines that the most probable root cause is the controller (94.2 %). It determines that the CPU (31.4 %), the VNF (31.4 %) or its configuration (31.4 %) could be the source of the failure. In case 2 (Fig. 7(b)), the actual root causes are three simultaneous link failures in the control link CL 1 and access links AL 1 and AL 2 .

The BN engine pinpoints those affected links CL 1 , AL 1 , AL 2 as the most probable root causes (31.1 %), having discarded the rest of elements. In the second scenario, we consider degradations in the SDN infrastructure affecting a service between H 1 and H 2 . This degradation may be explained by monitoring the CPU use on the network nodes, so the BN engine incorporates these observations. As consequence, the calculated root changes as a result of changes on these observations. We run the BN engine with two different CPU conditions on the nodes, shown in Fig. 8: (a) a non-loaded controller (CPU use 5%) with host H 2 heavily-loaded (CPU use 95%) and the rest of nodes with a CPU use between 5% and 95%), and (b) a heavily-loaded controller (CPU use 95%) and the rest of nodes with a CPU use between 5% and 95%. In situation (a), the BN engine determines that host H 2 is the most probable root cause (96.6%) due to its high CPU use, and so it discards all the links (1 %) as well as the controller (7.9 %) as probable causes. In the heavy load of CPU use in the controller (b), the BN engine selects the controller as the most probable root cause to explain the degradation on the SDN infrastructure (a transition of root cause probability from 7.9 % to 96.9 %).

Fig. 1 .

 1 Fig. 1. SDN-NFV scenario: (a) end-to-end service, and (b) network topology

Fig. 2 .

 2 Fig. 2. Self-Modeling based diagnosis framework

Fig. 3 .

 3 Fig. 3. Templates and Dependency Graphs: (a) network node and (b) link Each type of network node discovered in the network topology(controllers (C 1 … C n ), slave switches (SS 1 … SS n ), master switches (MS 1 … MS n ), hosts (H 1 … H n )) is an instance of this template.2) Link template The template of a network link (TL n Fig.3 (b)) is simpler and it is composed of the physical layer and one vertex. Each type of link (control links (CL 1 … CL n ), access links (AL 1 … AL n ), and inter switch links (IL 1 … IL n )) is an instance of this template.

Fig. 5 .

 5 Fig. 5. Dependency graphs of linear topology in: (a) out-of-band control, (b) in-band

Fig. 6 .

 6 Fig. 6. Speed as a function of the number of elements

Fig. 7 .

 7 Fig. 7. Root Cause analysis: (a) Case 1, (b) Case 2

Fig

  Fig. 8. Root Cause Analysis on changing CPU conditions

  The edge addition algorithm This algorithm adds the dependencies between the previously appended dependency graphs (GN n and GL n ) by the template instantiation algorithm. It receives as input the link descriptor and the global dependency graph (topologically sorted).  It puts one dependency graph of link GLn in between the two dependency graphs of the nodes GNn to be connected (these nodes are given in the link descriptor).  It then adds two E INTER edges from the dependency graph of the link GL n to the two dependency graph of nodes GN n . E INTER depicts the dependencies between GL n and GN n . In Fig.3(a), the dependency graph of one node GN n has two incoming E INTER edges from the dependency graphs of links GL i and GL j . In Fig.3(b), the dependency graph of one link GL n has two outgoing edges E INTER towards the dependency nodes of GN i and GN j .

	Algorithm: Topological Sorting
	IN: Global Dependency Graph G(VUNSORTED,EINTRA)
	OUT: Global Dependency Graph G(VSORTED,EINTRA)
	for each dependency graph appended to G
	for each layer in template
	obtain vertices of appended graph at current layer
	sort its vertices topologically
	end for
	end for
	III) Algorithm: Edge addition
	IN: Link Descriptor, Global Dependency Graph G(VSORTED,EINSIDE)
	OUT: Global Dependency Graph G(VSORTED,EINSIDE,EINTER)
	for each link in Link Descriptor
	extract end points attached to link
	for each end point in link
	EINTERadd edge from GLn[link,link] to GNn[node,card]
	end for
	end for

TABLE I

 I

	.	NUMBER OF VERTICES (V) AS A FUNCTION OF THE NUMBER
				OF HOSTS (NH)
	Topology/	4	8	16	32	64	128	256
	Tree	62	130	266	538	1082	2170	4346
	Linear	72	140	276	548	1036	2180	4356

  . 8. Root Cause Analysis on changing CPU conditionsVI. CONCLUSIONS AND FUTURE WORKThis paper considers solving two major problems towards self-diagnosis and resilient networks in the context of SDN and NFV: in such context it is needed to define a template or a model that describes the managed elements including physical, virtual infrastructure and other inner details such network cards, or CPU. To fill this gap, we define a template with finer granularity describing the essential managed elements within SDN and NFV. Furthermore, we specify and validate a self-modeling diagnosis that builds automatically at runtime the diagnosis model (dependency graph), which answers the challenges of updating the diagnosis model to identify the root causes. Our approach is suitable to any network topology and to any control type in SDN. In addition, it is independent from the controller implementation (e.g. Floodlight or OpenDaylight). Regarding future work, we will focus on the following points: -Reducing the uncertainty of the diagnosis: As a result of the finer granularity level of our proposed model, the uncertainty of the diagnosis is high. We foresee to adapt the methodology of Hounkonnou et. al.[START_REF] Hounkonnou | Active Self-Diangosis in Telecommunication Networks[END_REF] to reduce the uncertainty by updating the model progressively with observations obtained from new clients.-Learning network element templates to diagnose new root causes: we target to add learning mechanism to our selfmodeling based diagnosis. Statistical tests can be used as discussed in[START_REF] Bennacer | Optimization of fault diagnosis based on the combination of Bayesian Networks and Case-Based Reasoning[END_REF] to learn automatically node templates allowing the diagnosis of unexpected new root causes.-Modeling VNFs dependencies: Extension of the selfmodeling algorithm to model the VNF forwarding graphs that compose the service of each client online. I.

	Root Cause
	Probability (%):

ACKNOWLEDGMENTS

This work is partly funded by the French ANR under the ANR-14-CE28-0019 REFLEXION project


Amine Benelallam
email: amine.benelallam@inria.fr

Abel Gómez
email: abel.gomez@unizar.es

Massimo Tisi
email: massimo.tisi@inria.fr

ATL-MR: Model Transformation on MapReduce

Keywords: Model Transformation, Distributed Computing, Tool, ATL, MapReduce

The Model-Driven Engineering (MDE) paradigm has been successfully embraced for manufacturing maintainable software in several domains while decreasing costs and efforts. One of its principal concepts is rule-based Model Transformation (MT) that enables an automated processing of models for different intentions. The user-friendly syntax of MT languages is designed for allowing users to specify and execute these operations in an effortless manner. Existing MT engines, however, are incapable of accomplishing transformation operations in an acceptable time while facing complex transformations. Worse, against large amount of data, these tools crash throwing an out of memory exception.

In this paper, we introduce ATL-MR, a tool to automatically distribute the execution of model transformations written in a popular MT language, ATL, on top of a well-known distributed programming model, MapReduce. We briefly present an overview of our approach, we describe the changes with respect to the standard ATL transformation engine, finally, we experimentally show the scalability of this solution.

Introduction

Model-Driven Engineering (MDE) is gaining ground in industrial environments, thanks to its promise of lowering software development and maintenance effort [START_REF] Kärnä | Evaluating the use of domain-specific modeling in practice[END_REF]. It has been adopted with success in producing software for several domains such as the modernization of legacy software systems [START_REF] Bruneliére | MoDisco: A Model Driven Reverse Engineering Framework[END_REF]. Core concepts of MDE are the centrality of (software, data and system) models in all phases of software engineering and the au-tomation of model processing during the software life-cycle. Model Transformation (MT) languages have been designed to help users specifying and executing these model-graph manipulations. The AtlanMod Transformation Language (ATL) [START_REF] Jouault | ATL: A Model Transformation Tool[END_REF] is one of the most popular examples among them, and a plethora of transformations exist addressing different model types and intentions 1 .

In recent years, MDE has been witnessing the increasing complexity of systems and data that comes in MDE in the form of Very Large Models (VLMs) [START_REF] Clasen | Transforming Very Large Models in the Cloud: a Research Roadmap[END_REF]. Existing MDE tools, especially MT engines, are based on graph matching and traversing techniques, that makes them exposed to serious scalability issues in terms of memory occupancy and execution time. This stands in particular when MT execution is limited by the resources of a single machine.

With the broad availability of distributed clusters and programming models such as MapReduce [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF], a natural way to overcome these issues would be exploiting distributed systems for parallelizing model processing operations.

In this paper we introduce ATL-MR2 , a tool enabling the automatic execution of model transformations in ATL on top of MapReduce. The distribution is implicit, i.e. the syntax of the ATL is not modified and no primitive for distribution is added. Hence, no familiarity with distributed programming is required.

The paper is structured as follows. Section 2 introduces the syntax of ATL and its execution semantics by means of a running example. Section 3 describes an overview of our approach and details the implementation of our prototype engine. Section 4 discusses the evaluation results of our solution. Finally, Section 5 wraps up the conclusions and future works.

Background

In this section we give a brief introduction to ATL by means of model-driven reverse engineering case study. This case study was proposed as a simple but computationally expensive benchmark for MT engines [START_REF] Horn | The TTC 2013 Flowgraphs Case[END_REF]. We focus on one Fig. 1 shows an example of models, derived from a small program calculating a number factorial. Instructions are represented by rectangles, and variables by squares. An instruction points to the set of variables it defines or writes (def), and a set of variables it reads (use). The links cfNext and dfNext refer to the next control flow and data flow instructions respectively 3 . As it can be seen in the figure, the transformation changes the topology of the model graph.

The AtlanMod Transformation Language

In this paper we refer to an ATL implementation of the transformation named ControlFlow2DataFlow and available at the article website 2 . Listing 1 illustrates a rule of the transformation.

Model transformations in ATL [START_REF] Jouault | ATL: A Model Transformation Tool[END_REF] are unidirectional. They are applied to read-only source models and produce writeonly target models. Declarative rules abstract the relationship between source and target elements while hiding the semantics dealing with rule triggering, ordering, traceability management and so on.

ATL matched rules are composed of a source pattern and a target pattern. Both of source and target patterns might contain one or many pattern elements. Input patterns are fired automatically when an instance of the source pattern (a match) is identified, and produce an instance of the corresponding target pattern in the output model. Implicitly, transient tracing information is built to associate input elements to their correspondences in the target model.

Source patterns are defined as OCL [START_REF]Object Constraint Language, OCL[END_REF] guards over a set of typed elements, i.e. only input elements satisfying that guard are matched. In ATL, a source pattern lays within the body of the clause 'from'. For instance, in the rule SimpleStmt, the source pattern (Listing 1, lines 4-6) matches an element of type SimpleStmt that defines or uses at least a variable. Output patterns, delimited by the clause 'to' (lines 7-11) describe how to compute the model elements to produce when the rule is fired, starting from the values of the matched elements. E.g., the SimpleStmt rule produces a single 3 A detailed introduction can be found in the Dragonbook [START_REF] Lam | Compilers: Principles, Techniques, and Tools[END_REF], Chapter 9

Listing 1: ControlFlow2DataFlow -ATL rules (excerpt) ATL matched rules are executed in two phases, a match phase and an apply phase. In the first phase, the rules are applied over source models' elements satisfying their guards. Each single match, corresponds to the creation of an explicit traceability link. This link connects three items: the rule that triggered the application, the match, and the newly created output elements. At this stage, only the output pattern elements type is considered, bindings evaluation is left to the apply phase.

1 module ControlFlow2DataFlow ;
The apply phase deals with the initialization of output elements' features. Every feature is associated to a binding in an output pattern element of a given rule application. Features initialization is performed in two steps: 1) first the corresponding binding expression is computed, resulting in a collection of elements of the input model; 2) it is then passed to a resolution algorithm that, if needed, navigates the trace links produced in the match phase to determine the corresponding target elements.

MapReduce

MapReduce is a programming model and software framework developed at Google in 2004 [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF]. It allows easy and transparent distributed processing of big data sets while concealing the complex distribution details a developer might cross. Both Map and Reduce invocations are distributed across cluster nodes, thanks to the Master that orchestrates jobs assignment. Input data is partitioned into a set of chunks called Splits. Every split comprises a set of logical Records.

Given the number of Splits and idle nodes, the Master node decides the number of workers (slave machines) for the assignment of Map jobs. Each Map worker reads one or many Splits, iterates over the Records, processes the records and stores the results locally. When Map workers finish, the Master forwards these locations to the Reduce workers that sort the records by key so that all occurrences of the same key

ATL-MR

Our ATL-MR transformation engine follows a data-distribution scheme, where each one of the nodes that are computing in parallel takes charge of transforming a part of the input model. Its source code is available at the paper's website.

The engine is built on top of the ATL Virtual Machine (EMFTVM [START_REF]ATL EMFTVM[END_REF]) and Apache Hadoop [START_REF]Apache Hadoop[END_REF]. Fig. 2 shows how the ATL transformation of our running example is executed on top of three nodes, two map and one reduce workers. The input model is equally split over map workers. In the map phase, each worker runs independently the full transformation code but applies it only to the assigned subset of the input model. We call this phase Local matchapply. Afterwards each map worker communicates the set of model elements it created to the worker that performs the reduce phase, together with trace information. These trace links (grey arrows in Fig. 2) encode the additional information that will be needed to resolve the binding, i.e. identify the exact target element that has to be referenced based on the tracing information. The reduce worker is responsible of gathering partial models and trace links from the map workers, and updating properties value of unresolved bindings. We call this phase Global resolve. Each node in the system executes its own instance of the ATL VM but performs either only the local match-apply or the global resolve phase.

At the map phase, input splits are assigned to map workers. Each one of these splits contains a subset of the input model for processing. However, each worker has a full view of the input models in case it needs additional data for bindings computation. For example, executing the binding dfNext over the method fact(int a), results in {while(a)>0, r*=a--;} (grey arrows in Fig. 2). Since while(a) and fact(int a) reside in the same node, a dfNext reference between them is created in the target model during this map phase. At the beginning of the reduce phase, all the target elements are created, the local bindings are populated, and the unresolved bindings are referring to the source elements to be resolved. This information is kept consistent in the tracing information formerly computed and communicated by the mappers. Then it resolves the remaining reference bindings by iterating over Figure 3: ATL-MR's integration within MONDO platform the trace links. For each trace link, the reducer iterates over the unresolved elements of its property traces, resolves their corresponding element in the output model, and updates the target element with the final references. In the right-hand side of Fig. 2 all the trace properties have been substituted with final references to the output model elements.

The standalone and distributed ATL engines share most of the code. In the standard ATL VM, the transformation engine iterates over the matched rules, and looks for the elements that match its application condition (guard). Instead, our distributed VM iterates over each input model element, and checks if it is matched by an existing rule. In this perspective, we extended the ATL VM with a minimal set of functions allowing the VM to run either in standalone or distributed mode. In particular, the distributed VM is required to factorize and expose methods for launching independently small parts of the execution algorithms, for instance the transformation of single model elements. More information about the tool usage and deployment can be found at the tool's website 2 .

Platform Integration

Besides running independently, ATL-MR can be used as a component of MONDO, an integrated platform to support scalable MDE [START_REF][END_REF]. MONDO consists of a distributed platform and an Eclipse-based tool workbench counterpart through which developers can interact with the facilities provided by the platform. As depicted in Fig. 3 ATL-MR client is implemented as OSGI bundle integrated with the MONDO client. It comes with a command line-based API enabling four primary operations, namely launch and kill for launching and stopping a transformation respectively, whereas, status returns the status of a transformation job. Finally, listJobs, returns the list of running jobs. On the server side, ATL-MR is implemented as part of the MONDO server by means of a web servlet communicating with the Resource Manager of the Hadoop cluster.

Limitations

Currently, our ATL VM supports only the default EMF serialization format XMI (XML Metadata Interchange). This file-based representation faces many issues related to scalability. In particular, models stored in XMI need to be fully loaded in memory, but more importantly, XMI does not support concurrent read/write. This hampers our tool at two levels, first, all the nodes should load the whole model even though if they only need a subset of it. This prevents us from transforming very big models that would not fit in memory. The second one concerns the reduce phase parallelization, and this is due to the fact that only one mapper can write to the output XMI file at once. In a recent work, we extended an existing persistence backend NeoEMF [START_REF] Gómez | Map-based transparent persistence for very large models[END_REF] with support for concurrent read/write and on-demand loading [START_REF] Gómez | Decentralized Model Persistence for Distributed Computing[END_REF] on top of Apache HBase [START_REF]Apache HBase[END_REF]. In our current work, we are coupling it with our VM to solve these two particular issues.

Experimentation

In this section we experimentally evaluate the scalability of our approach by conducting two different but complementary experiments. We run our experiments in our running example and compare how it performs in two different test environments (clusters). The transformation covers a sufficient set of declarative ATL constructs enabling the specification of a large group of MTs. We use as input different sets of models of various sizes, reverse-engineered from a set of automatically generated Java programs. While the first one shows a quasi-linear speed-up w.r.t. the cluster size for input models with similar size, the second one illustrates that the speed-up grows with increasing model size.

For the first experiment we have used a set of 5 automatically generated Java programs with random structure but similar size and complexity. The source Java files range from 1 442 to 1 533 lines of code and the execution time of their sequential transformation ranges from 620s to 778s. The experiments were run on a set of identical Elastic MapReduce clusters provided by Amazon Web Services. All the clusters were composed by 10 EC2 instances of type m1.large. Each execution of the transformation was launched in one of those clusters with a fixed number of nodes -from 1 to 8 -depend- To investigate the correlation between model size and speed-up we execute the transformation over 4 artificially generated Java programs with identical structure but different size (from 13 500 to 105 000 lines of code). This time the experiments have been executed in a virtual cluster composed by 12 instances built on top of OpenVZ containers running Hadoop 2.5.1.

As shown in Fig. 5, the curves produced by Experiment II are consistent to the results obtained from Experiment I, despite the different model sizes and cluster architectures. Moreover, as expected, larger models produce higher speedups: for longer transformations the parallelization benefits of longer map tasks overtakes the overhead of the MapReduce framework.

Conclusion and Future Work

In this paper we introduce ATL-MR a tool to enable the distribution of ATL transformation in a distributed environment. We show its interoperability by integrating it with a distributed MDE platform, MONDO. Also, we experimentally show the good scalability of our solution. Thanks to our publicly available execution engine, users may exploit the availability of MapReduce clusters on the Cloud to run model transformations in a scalable and fault-tolerant way.

In our future work we plan to improve the efficiency of our tool, by addressing related research aspects. We aim to investigate: (i) coupling with the transformation engine a distributed model-persistence backend supporting concurrent read/write, and on-demand loading, (ii) efficiently distributing the input model over map workers with the aim to optimize load balancing and minimize workload.

Figure 1 :

 1 Figure 1: ControlFlow2DataFlow transformation example

2 c 3 r u l e SimpleStmt { 4 from 5 s 7 t o 8 t 11 ) 12 }

 2345781112 r e a t e OUT : DataFlow from IN : ControlFlow ; : ControlFlow ! SimpleStmt ( not ( s . def -> 6 isEmpty ( ) and s . use->isEmpty ( ) ) ) : DataFlow ! SimpleStmt ( 9 txt <s . txt , 10 dfNext <s . computeNextDataFlows ( ) element of type SimpleStmt. A set of OCL bindings specify how to fill each of the features (attributes and references) of the produced elements. The binding at line 9 copies the textual representation of the instruction, the binding at line 10 fills the dfNext link with values computed by the computeNextDataFlows OCL helper.

Figure 2 :

 2 Figure 2: ControlFlow2DataFlow example on MapReduce

Figure 4 :

 4 Figure 4: Speed-up obtained in Experiment I

Figure 5 :

 5 Figure 5: Execution times and speed-up on Experiment II

http://www.eclipse.org/atl/atlTransformations/

https://github.com/atlanmod/ATL_MR/

Acknowledgments

This work is partially supported by the MONDO (EU ICT-611125) project. Many thanks to Antonio García-Domínguez for his valuable help while integrating ATL-MR with the MONDO platform.