G Legrain

N Moës Gem

N Moës

Adaptive anisotropic integration scheme for high-order fictitious domain methods: Application to thin structures

Keywords: X-FEM, high-order, fictitious domain, anisotropic meshing, solid-shell, shell, plates, beams

la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fictitious domain methods have recently gained more attention because of their favourable properties to conduct design-through analysis [START_REF] Schillinger | An isogeometric design-through-analysis methodology based on adaptive hierarchical refinement of NURBS, immersed boundary methods, and Tspline CAD surfaces[END_REF]. Indeed, this concept aims at creating more efficient part designs thanks to the tight relationship between CAD and analysis. It allows designers to draw parts with mechanical efficiency in mind. Such strategies require a seamless transition between geometrical representation (CAD surfaces for instance) and analysis. Isogeometric analysis [START_REF] Cottrell | Isogeometric Analysis[END_REF] was proposed with this objective in mind, and has proved to be a very efficient method. However, meshing is still necessary: it requires clean CAD surfaces, a volumetric parametrization of the solid [START_REF] Cottrell | Isogeometric Analysis[END_REF] (whereas the CAD only deals with surface representation), and also trimming-free surfaces [START_REF] Kim | Isogeometric analysis for trimmed CAD surfaces[END_REF][START_REF] Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF] , although this last topic is the subject of numerous recent works [START_REF] Beer | A simple approach to the numerical simulation with trimmed CAD surfaces[END_REF][START_REF] Kim | Isogeometric analysis for trimmed CAD surfaces[END_REF][START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Kollmannsberger | Parameter-free, weak imposition of Dirichlet boundary conditions and coupling of trimmed and non-conforming patches[END_REF][START_REF] Nagy | On the numerical integration of trimmed isogeometric elements[END_REF]. Thanks to their ability to use unfitted meshes, fictitious domain methods are also an answer to this objective. Albeit known for decades [START_REF] Saulèv | On the solution of some boundary value problems on high performance computers by fictitious domain method[END_REF], fictitious domain methods have been recently revisited in the context of the X-FEM [START_REF] Belytschko | Structured extended finite element methods for solids defined by implicit surfaces[END_REF][START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: Application to homogenization[END_REF][START_REF] Lian | Image-based computational homogenization and localization: comparison between X-FEM/levelset and voxel-based approaches[END_REF][START_REF] Legrain | Treatment of nearly-singular problems with the X-FEM[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF][START_REF] Legrain | Numerical simulation of CAD thin structures using the eXtended Finite Element Method and Level Sets[END_REF], Cartesian grids [START_REF] Nadal | Efficient Finite Element Methodology Based on Cartesian Grids: Application to Structural Shape Optimization[END_REF] for low-order approximations, and also in the context of higher order approximations with the introduction of the Finite Cell method [START_REF] Düster | The finite cell method for threedimensional problems of solid mechanics[END_REF][START_REF] Rank | Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures[END_REF][START_REF] Parvizian | Finite cell method -h and p extension for embedded domain problems in solid mechanics[END_REF][START_REF] Düster | Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method[END_REF][START_REF] Schillinger | Small and large deformation analysis with the p-and B-spline versions of the Finite Cell Method[END_REF][START_REF] Schillinger | An unfitted hp-adaptive finite element method based on hierarchical B-splines for interface problems of complex geometry[END_REF]. This method can be seen as a hybrid between p-fem and fictitious domains, retaining advantages of both (i.e. exponential convergence for smooth solutions without meshing burden). The difficulty shared by these methods lies in the cost and accuracy of the integration of the weak form, as meshing difficulties are transferred into the proper integration of this quantity. For low order approximations, the geometrical error is usually of higher order than the approximation error: "coarse" polygonal discretizations are sufficient if their size is similar to the size of the approximation elements. Unfortunately, it is not the case for higher-order approximations where the geometrical error has to be kept smaller than the approximation error (which decreases exponentially). Sub-grid integration strategies have been proposed in the context of both Finite Cell [START_REF] Parvizian | Finite cell method -h and p extension for embedded domain problems in solid mechanics[END_REF] and higher-order X-FEM [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF]. In the first case, integration points in the sub-grid are filtered based on their location whereas a fine triangulation of the surface is built for the latter. These two approaches are satisfactory, but can require quite large computational resources for the integration due to the large number of integration points, unless pre-computed strategies like in [START_REF] Yang | An efficient integration technique for the voxel-based finite cell method[END_REF] are considered. Recent contributions [START_REF] Mousavi | Numerical integration of polynomials and discontinuous functions on irregular convex polygons and polyhedrons[END_REF][START_REF] Sudhakar | Quadrature schemes for arbitrary convex/concave volumes and integration of weak form in enriched partition of unity methods[END_REF][START_REF] Müller | Highly accurate surface and volume integration on implicit domains by means of moment-fitting[END_REF][START_REF] Joulaian | Numerical integration of discontinuities on arbitrary domains based on moment fitting[END_REF] propose to use moment fitting techniques [START_REF] Lyness | Moderate Degree Symmetric Quadrature Rules for the Triangle[END_REF][START_REF] Lyness | Quadrature Rules for Regions Having Regular Hexagonal Symmetry[END_REF][START_REF] Dunavant | High degree efficient symmetrical Gaussian quadrature rules for the triangle[END_REF] to construct integration rules on the fly. Blended integration cell methods have also been proposed in order to decrease the integration cost thanks to large mapped integration cells [START_REF] Strouboulis | The Generalized Finite Element Method: An Example of its Implementation and Illustration of its Performance[END_REF][START_REF] Legrain | A NURBS enhanced extended finite element approach for unfitted CAD analysis[END_REF][START_REF] Marco | Exact 3D boundary representation in finite element analysis based on Cartesian grids independent of the geometry[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary methods[END_REF][START_REF] Kudela | Smart octrees: accurately integrating discontinuous functions in 3D[END_REF][START_REF] Fries | Higher-order accurate integration of implicit geometries[END_REF]. All these methods can improve the efficiency of the integration although the most robust method remains the sub-grid based integration. In this paper, we propose an efficient variant of sub-grid integration rules that improves its numerical efficiency in low-curvature areas, and also evaluate new possibilities opened by this method. In particular, we will focus on the treatment of slender structures that are typically challenging for this family of methods, as the size of the integration cells must be smaller than the thickness of the structure. The method is based on an anisotropic mesh adaptation strategy [START_REF] Frey | Mesh Generation: Application to Finite Elements[END_REF] which consists in optimizing not only mesh density, but also elements shape. A large literature is devoted to this topic, interested reader is referred to [START_REF] Huang | Anisotropic Mesh Adaptation for Variational Problems Using Error Estimation Based on Hierarchical Bases[END_REF][START_REF] Piggott | Anisotropic mesh adaptivity for multi-scale ocean modelling[END_REF][START_REF] Picasso | Numerical study of the effectivity index for an anisotropic error indicator based on Zienkiewicz-Zhu error estimator[END_REF][START_REF] Bui | An accurate anisotropic adaptation method for solving the level set advection equation[END_REF][START_REF] Bourgault | Anisotropic error estimates and space adaptivity for a semidiscrete finite element approximation of the transient transport equation[END_REF][START_REF] Farrell | An anisotropic Zienkiewicz-Zhu-type error estimator for 3D applications[END_REF][START_REF] Ducrot | Anisotropic Level Set Adaptation for Accurate Interface Capturing[END_REF][START_REF] Abgrall | An immersed boundary method using unstructured anisotropic mesh adaptation combined with level-sets and penalization techniques[END_REF][START_REF] Quan | Anisotropic mesh adaptation with optimal convergence for finite elements using embedded geometries[END_REF] to name a few. The adaptation strategy that is considered here is relying on the metric tensor proposed by Coupez [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Coupez | Implicit Boundary and Adaptive Anisotropic Meshing[END_REF]. This approach is based on the level-set representation of the interface and ensures a proper representation of this surface with a minimum number of elements. The resulting mesh is naturally refined near large curvature regions and elongated along low curvature zones. The proposed method enables to decrease the number of integration points, but can also effectively capture the very thin areas of the domain.

We thus propose to study the performances of fictitious domain methods for this family of problems, which are generally approximated by means of shell and plate elements. These models are computationally efficient as they require very few dofs through the discretization of the mid-surfaces. However this meshing step can be tedious and these elements are unable to represent complex through the thickness behaviours, and lead to modelization issues for contact and in plate to solid transition zones. Thereby leading researchers to develop so-called solid-shell elements that are based on displacement-only kinematic, with one layer of element across the thickness in order to improve their computational efficiency. Such methods can effectively model thick and thin areas, large strain problems, and capture through the thickness phenomena (especially for nonlinear problems). The meshing step is also alleviated (although not removed). Naive solid-shell elements are unfortunately subjected to membrane and shear locking so that adapted formulations relying generally on assumed strain formulations with reduced integration and hourglass control are usually used [START_REF] Abed-Meraim | SHB8PS--a new adaptative, assumed-strain continuum mechanics shell element for impact analysis[END_REF][START_REF] Alves De Sousa | A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness: Part I-geometrically linear applications[END_REF][START_REF] Reese | A large deformation solid-shell concept based on reduced integration with hourglass stabilization[END_REF][START_REF] Bassa | A new nine-node solid-shell finite element using complete 3D constitutive laws[END_REF]. This family of methods have been recently extended in the context of isogeometric analysis by Bouclier et al. [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF][START_REF] Bouclier | Efficient isogeometric NURBS-based solid-shell elements: Mixed formulation and Bbar-method[END_REF][START_REF] Bouclier | An isogeometric locking-free NURBSbased solid-shell element for geometrically nonlinear analysis[END_REF]. In addition, note also the contribution of Düster et al. [START_REF] Düster | pq-Adaptive solid finite elements for threedimensional plates and shells[END_REF] who proposed effective p-fem formulations for plates en shells. All these methods rely on the meshing of the structure which can be tedious for complex industrial parts. We therefore propose to evaluate the performances of high-order fictitious domain methods in this context in order to combine both the performances of fictitious domains (meshing ease) and high-order methods (exponential convergence). The proposed method can be seen as an extension of the method proposed in [START_REF] Legrain | Numerical simulation of CAD thin structures using the eXtended Finite Element Method and Level Sets[END_REF].

The remainder of the paper is organized as follows: first, the anisotropic mesh adaptation method used in this contribution is introduced and illustrated. Then, its coupling with fictitious domain methods is presented. In a second part, the strategy is assessed both with respect to the integration cost and with respect to solid-shell applications. The method is finally verified against common shell benchmarks in both 2D and 3D.

2 Fictitious domain strategies based on anisotropic mesh adaptation

Motivation

The use of fictitious domain methods for thin structures can be traced back to only two contributions [START_REF] Rank | Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures[END_REF][START_REF] Legrain | Numerical simulation of CAD thin structures using the eXtended Finite Element Method and Level Sets[END_REF]. The Shell Finite Cell [START_REF] Rank | Shell Finite Cell Method: A high order fictitious domain approach for thin-walled structures[END_REF] consists in mapping a 2D parametric space onto a thin 3D space. The Finite-Cell method is used in the parametric space to streamline the representation of complex geometrical features. In [START_REF] Legrain | Numerical simulation of CAD thin structures using the eXtended Finite Element Method and Level Sets[END_REF], a moderate order (quadratic) X-FEM approximation in the physical space was used. This approach requires a large number of approximation elements in order to be able to resolve the thickness of the structures and represent geometrical edges with a good accuracy (see Figure 1 (a)). Increasing the order of the approximation and using sub-grid level-sets [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF] was also advocated, but led to a large amount of integration points (see Figure 1 (b)). Indeed, the size of the geometrical cells must be smaller than the thickness of the structure so that it can be represented. In the case of mixed slender and massive zones, it either leads to massive integration costs or the disappearance of the slender zones, depending on the user-defined size of these cells. In addition, the edges of the geometry are lost as their scale is below the geometrical resolution of the geometrical mesh.

The objective here is to propose an alternative geometrical representation strategy for such applications. This approach, which is based on anisotropic mesh adaptation, is top to bottom and allows to capture both small and large geometrical features. It should also be very efficient for representing low-curvature areas that are common for industrial designs.

Figure 1: Strategies for thin structures (thick lines: approximation elements, dashed lines: geometrical elements): (a) moderate order X-FEM [START_REF] Legrain | Numerical simulation of CAD thin structures using the eXtended Finite Element Method and Level Sets[END_REF]; (b) sub-grid level-set [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF][START_REF] Legrain | High order X-FEM and levelsets for complex microstructures: Uncoupling geometry and approximation[END_REF].

Anisotropic mesh adaptation for level-sets

The objective of anisotropic mesh adaptation is to construct a mesh for which the approximation error is uniform and below a given threshold. The main difference with classical (isotropic) mesh adaptation lies in the shape of the elements. The only parameter for isotropic mesh refinement is the element size field across the structure of interest, whereas anisotropic mesh refinement allows to specify the size of the elements in all the directions to achieve error uniformity. If the solution has a slow evolution in one direction, then elements size can be increased along this direction without any particular influence on the error level. In the example presented in Figure 2, the error is obviously lower for the anisotropic mesh for a given number of vertices.

The method relies on the construction of a metric tensor, M(x) which is used by the mesher to redefine the way lengths are computed depending on the direction. This metric tensor is usually obtained on an element-basis [START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF], but in practice it is more convenient to work with nodal metrics as it leads to continuous metric fields in space. These continuous fields greatly simplify metric interpolation and extrapolation operations. The approach considered in this contribution follows this alternative path, and relies directly on a nodal metric (without any post-processing of elementary quantities). This approach was proposed in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF], and depends on a so called length-distribution tensor which is derived from a statistical concept introduced in polymer science [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF]. The originality lies in the construction of the metric, and the associated edge-based error estimator which avoids complex Hessian reconstructions like those proposed in [START_REF] Abgrall | An immersed boundary method using unstructured anisotropic mesh adaptation combined with level-sets and penalization techniques[END_REF][START_REF] Claisse | Levelsets and anisotropic mesh adaptation[END_REF]. In the case of level-sets, the target region for the adaptation process is not the whole domain, but the zone near the iso-zero. This is achieved here by working with the regularized Heaviside function of the level-set: This function is a proper choice as it is C ∞ and concentrates the interpolation error solely in the vicinity of the iso-zero. The features of the metric are now illustrated for some 2D examples. We only focus on some qualitative illustrations, in order to highlight the main features of the adaptation procedure, as it has already been deeply studied in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Coupez | Solution of high-Reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF][START_REF] Coupez | Implicit Boundary and Adaptive Anisotropic Meshing[END_REF]. To this end, different geometries are considered: (i) a plane interface, (ii) a circular interface and (iii) a corner interface. The first geometry involves no curvature, the second has a constant one, and the third one a varying curvature. Adapted meshes for different values of regularization parameter ε are presented in Figure 3.

u ε (x) = ε tanh φ(x) ε . (1)
As expected, the elements are aggregated in a band of length 2 √ 2ε and elongated along the directions with small curvature. In regions of large curvature (near the apex of the corner), then the elements are isotropic. In addition, it should be noted that decreasing ε allows to master the geometrical accuracy of the discrete geometry, as the approximated iso-zero lies in the 2 √ 2ε band around the real iso-zero.

Anisotropic mesh adaptation for fictitious domains

The anisotropic integration mesh which is generated by the anisotropic mesh adaptation procedure has to be linked to the approximation mesh, so that the weak formulation can be integrated properly. In the case of sub-grid level-sets or with the Finite-Cell, this connection step is trivial as the integration cells are embedded within the approximation mesh (see Figure 1 (b)). In the present case, an ad-hoc strategy has to be proposed as the anisotropic integration mesh is completely unrelated to the approximation mesh (see Figure 4(a)).

In order to build a set of embedded integration cells from the anisotropic mesh, we choose to focus on regular quadrangular or hexahedral approximation meshes. Then, the regular pattern of this approximation mesh is used to slice the integration mesh along the boundaries of the approximation elements, thereby leading to a set of embedded anisotropic integration cells ((see Figure 4(b)). This operation can be done robustly and efficiently by re-using level-set slicing routines. Note however that this strategy is not restricted to regular grids and could be easily extended to octree-based approximation meshes, applying it recursively in the leaves of the octree. In a last step (Figure 4(c)), integration cells are only kept for the approximation elements containing the interface to be captured. This allows to restrict integration effort only where it is needed. After this procedure, one obtains a suitable integration mesh for the fictitious domain procedure.

ε = 10 -1 ε = 10 -2 ε = 10 -3 ε = 10 -3 (zoom)
Note that this slicing procedure leads to an increase in the number of integration cells, which could be balanced by only retaining the narrow-band of the anisotropic integration mesh. The development of such a strategy is currently in progress. Throughout this contribution, all the integration cells are equipped with full Gauss quadrature rules, designed so that the polynomial quantities involved in the formulation can be evaluated exactly. This leads to a large number of integration points, but one could also use a non-uniform quadrature rule, such as the one proposed in [START_REF] Abedian | Performance of Different Integration Schemes in Facing Discontinuities in the Finite Cell Method[END_REF]. It consists in lowering the number of integration points for the very small cells that are created near the interfaces in order to further improve the numerical efficiency of the integration.

Comparison with octree adaptation

This section is devoted to the characterization of the strengths and weaknesses of the anisotropic procedure compared to the classical octree-based strategy, in terms of geometrical accuracy only. The influence of this geometrical accuracy on the solution of the physical problem is studied in Section 3.1. The two methods are compared for fixed level-sets, by measuring their geometrical accuracy with respect to the number of nodes in the meshes they produce. Given a level-set φ(x) and a mesh, the (absolute) geometrical error E g associated to the mesh is defined as the integral of the distance between the discrete iso-zero of the level-set Γ h and the exact iso-zero. This quantity can be evaluated easily when the exact level-set of the problem is known:

E g = Γ h |φ(x)|dΓ. (2)
We study the compared evolution of this error between octree and anisotropically adapted meshes. Throughout this section, octree and anisotropic results are depicted with solid and dashed lines respectively. Moreover, these curves share similar ticks for a given geometry. First, circular interfaces are considered with increasing radii, from 0.05 to 5 (see Figure 5 (a)). Octree meshes were generated by changing their depths from 5 to 10, and the anisotropic meshes were constructed with the same target number of nodes. it is not influenced by the curvature of the curve. This can be explained by the fact that in this case, the error is already homogeneous around the interface for isotropic meshes (the anisotropic strategy also tends to create isotropic element in this case, see Figure 3). Then, the case of a geometry mixing low and large curvature is considered in Figure 6 (a). The interface is obtained by the intersection of two off-centred circular level-sets of varying radii. The corresponding interface tends to a straight corner when the radius of the circles is increased. From the convergence study presented in Figure 6 (b), we can draw the following conclusions: (i) The presence of large and low curvature areas it better handled by the anisotropic approach. This stems from the fact that the anisotropic adaptation method is able to optimize the mesh when the geometrical error is non uniformly distributed in the domain, whereas the uniform size of the elements generated by the octree cannot. low and large curvature areas, the better the performances of the anisotropic method.

(iii) The efficiency of the anisotropic algorithm is at least on par with the octree scheme.

Finally, note that the performances of the anisotropic strategy are hindered by the layers of elements in the zones near the interface, as they do not contribute to the improvement of the geometrical accuracy. Methods making use of the narrow-band only may further improve its performances (more than for the octree).

Numerical examples

The numerical examples are divided in two parts: in the first one, the performance of the proposed approach is assessed with respect to the ratio between the cost of the numerical integration and the accuracy. In a second part, the potential of the use of such approach is highlighted in the case of the treatment of thin structures.

Integration cost

In the previous section, the parameters that influence the geometrical performances of the anisotropic adaptation scheme were studied. We focus now on the implications of these conclusions when dealing with physical problems. The objective of this section is to study to what extent do the geometrical conclusions extend to full calculations.

Circular Interface

Consider an infinite plate with a circular hole in its center which is loaded in uniaxial tension σ ∞ = 1.0 MPa along x axis (see Figure 7). The analytical solution of this prob-lem is given in [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF] and is used here as a reference solution for error monitoring. Only a square of length L = 2.0 mm with a circular hole of radius a = 0.55 mm at its center is considered here. Finally, tractions computed from the exact solution are applied on the boundary of the domain, and three dofs are prescribed to zero in order to avoid rigid body motion. The exact solution of this finite problem is therefore the same as the exact solution inside the boundary. Young's modulus is set to 1.0 MPa, and Poisson's ratio to 0.3. A p-convergence study is performed for p = 1 to 12, and the error in the energy norm is monitored with respect to the number of dofs and integration points. The approximation mesh is composed of 2 × 2 quadrangular elements with Legendre hierarchical shape functions. Typical integration meshes are depicted in Figure 8. In the octree case, the influence of the quality of the integration is presented in Figure 9 with h g = h/16 and h g = h/64 near the interface (where h g stands for the size of the elements in the geometrical mesh crossed by the interface). This typical convergence curve highlights the interplay between geometrical and approximation errors. We now restrict ourselves to the coarsest mesh, and these results are compared with those obtained using anisotropic integration meshes. These meshes were designed so that the number of integration cells is similar between the two methods. With such coarse meshes (see Figure 10), the anisotropic regularization parameter is set to 3. 10 -2 . The convergence curves corresponding to these three integration meshes are depicted in Figure 11. If one compares the convergence with respect to the number of dofs, then one can see that the behaviour of the anisotropic integration mesh is similar to the one of the octree, even if better error levels can be obtained. More interesting is the convergence with respect to the number of integration points, as these curves represent the integration effort which is required to attain a given error level. It is shown that for a similar number of integration points, the anisotropic integration performs better than the octree one (see the curve "Aniso 259"): for the same number of integration points, the asymptotic error is divided by 2. Before the point where the geometrical error exceeds the approximation error, the two schemes perform the same. More interestingly, the coarsest anisotropic mesh (curve "Aniso 212") constantly outperforms the octree one: it is more efficient and gives more accurate results than the octree. These results, however, do not show the full potential of the method: Indeed, it was shown in Section 2.4 that in the case of a constant curvature, the performances of the anisotropic strategy are similar to the performances of the octree. Hence, a case with non-uniform curvature is considered in the next section.

Varying curvature interface

Consider the example presented in Figure 12. It represents a square of length 2 × 2 with a square hole of length a = 0.8 at the center of the frame. A Laplace equation is solved on this domain: source and flux are prescribed respectively in the domain and on the boundary so that the solution is u(x, y) = cos(πx/a) sin(πy/a) (see Figure 12 for a sketch). The problem is solved by means of a 4 × 4 mesh which is presented in Figure 12. The geometry is either represented using an octree mesh or using an anisotropic mesh.

A first solution is obtained by means of an octree mesh with h g = h/32: in this case, the geometry is represented exactly, as the corners of the square are located on the edges of the elements. A p-convergence is conducted, from p = 1 to 8, and it can be seen in Figure 14(a) that an exponential convergence is obtained in the energy norm. This is consistent, since the geometry is captured exactly and the solution is regular. The center of the hole is now moved to (-0.18, -0.09), so that it is not located on the edges of the elements.

Figure 13: Integration cells for an octree integration scheme with h g = h/32.

The corresponding integration mesh is shown in Figure 13 (1168 nodes): the error in the perimeter length of the square is 0.5% in this case. It can be seen in Figure 14(a) that this geometrical error causes a severe loss of accuracy as the convergence stagnates at a level of 0.4%.

Anisotropic integration is now considered (with a regularization parameter of 10 which leads to the 639 nodes mesh presented in Figure 15. Using this mesh for the geometry enables to recover the exponential convergence for all the polynomial range. However, we are also interested in the efficiency of the integration scheme. Various integration schemes are compared in Figure 14(b), where the error in the energy norm is plotted with respect to the number of integration points. One can see that if an octree mesh is used, then error levels similar to the anisotropic mesh can be obtained only by selecting very fine elements in the vicinity of the corners (h/256). The price to pay for the integration is an increase of one order of magnitude of the number of integration points. If the target error level is near 10 -3 , then coarser octree meshes can be selected (h g = h/64 is sufficient here). Nonetheless, the same error level can be obtained if an anisotropic integration mesh is used: in this case, the number of integration points can be reduced by one order of magnitude. Similar conclusions can be drawn for larger target error levels.

In conclusion, we have demonstrated that the use on the anisotropic integration strategy enables to save a large number of integration cells if the interface contains both low and large curvature zones (which is the case for a large number of industrial structures). On the contrary, the benefits of the anisotropic strategy are smaller when constant curvature zones are considered (see Sections 2.4 and 3.1.1). Yet, in all cases the strategy is at least as efficient as octree integration. In addition, the method works "top to bottom" in contrast to the octree: it enables to capture small features that could be lost by the octree (below the resolution of the maximal depth).

2D Thin structures

Thanks to its features, the proposed method is well adapted if thin structures are of interest. As these slender geometries can be represented accurately and efficiently, we try to assess fully three-dimensional fictitious domain strategies for such structures. The salient benefits are the ease of meshing, the possibility to handle true three-dimensional solutions with simple transition between thin and massive zones and general nonlinear models. As locking is the main difficulty for solid-shell methods, we focus on the robustness of the method with respect to locking by applying the method to classical test-cases from the literature.

Quarter circle beam

Consider a quarter circle beam as depicted in Figure 16(a). If the width of the beam is assumed to be small, then the problem can be solved in 2D (on the cross-section) with the plane stress assumption. Such geometry is prone to locking when discretized with classical solid finite elements. The beam is blocked at its bottom face and a unit vertical load is distributed on the left face. The radius of the beam is kept fixed whereas the thickness h is set to R/10, R/100 or R/1000, which correspond respectively to slenderness ratios of 10, 100 and 1000. Solving this problem with a Timoshenko beam model gives the following deflection:

u ref y = F R h 2 µ b + 3 F π R 3 E b , ( 3
)
where µ is the shear modulus, E Young's modulus and b the width of the beam. The problem is solved by means of fictitious domain method with the anisotropic integration algorithm. The resulting integration meshes for the slenderness ratios given above are presented in Figure 17. All meshes were designed with a regularization parameter ε = 10 -6 and a target of 1100 nodes for the two smaller thickness ratios, and 1500 nodes for the thinnest geometry. The approximation mesh is kept fixed with 16×16 hierarchical elements in the bounding box of the beam. The corresponding remaining cells are presented in Figure 18. A p convergence is conducted with p from 1 to 8, and the vertical displacement on the 19. It can be seen that the convergence of the displacement towards the reference solution is very fast, and almost not influenced by the slenderness ratio. Note also that such a problem could be solved by means of conforming high-order elements [START_REF] Düster | The p-version of the finite element method for three-dimensional curved thin walled structures[END_REF] or isogeometric elements [START_REF] Bouclier | On the development of NURBS-based isogeometric solid shell elements: 2D problems and preliminary extension to 3D[END_REF].

y x R h F A U A
Figure 19: Convergence of the displacement for slenderness ratios of 10, 100 and 1000.

A discussion on conditioning

The influence on the condition number of the finite element operator is now discussed. Consider the problem of a clamped beam of width t depicted in Figure 20a. The beam is immersed in a mesh composed of 16 × 1 or 16 × 3 quadrangular elements, depending on its slenderness. A hierarchical Legendre approximation is defined on these elements, and the evolution of the condition number of the operator against the increase of the slenderness of the beam L/t is monitored during p-refinement (p = 1 • • • 6). The condition number of the operator κ(A) = λ max /λ min is computed through the evaluation of the extreme eigenvalues of A, λ max and λ min . These eigenvalues are computed by means of an implicitly restarted Arnoldi iterations algorithm [START_REF] Arnoldi | The principle of minimized iterations in the solution of the matrix eigenvalue problem[END_REF]. Note that a Shift and Invert algorithm is used to estimate efficiently λ min : Rather than finding this minimal eigenvalue, the following modified eigenvalue problem is considered: (A -σI) -1 x = νx, where σ is close to λ min , and ν = 1/(λ min -σ). This strategy enhances the convergence of the algorithm to the desired portion of the spectrum. The evolution of the condition number is depicted in Figure 21 (a). The impact of the slenderness ratio L/t is clearly seen, as although large slenderness can be considered (up to 1000), the condition number increases at a rate which is clearly linked to L/t. For slenderness ratios larger 200, round-off errors prevent to compute accurately λ min when p > 4 (this is due to the inverse of A that appears in the Shift and Invert eigenvalue problem), so that the polynomial range for plotting is reduced. It has to be noted that p = 4 led to converged (in 2D) or almost converged results (in 3D) for all the numerical examples presented in this paper. It is also interesting to mention that slenderness is not the only cause of these conditioning issues: massive beams can also lead to very large condition number (see Figure 20b for the example of a beam with L/t = 7.9). This behaviour is discussed carefully in [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF] where the authors have shown that the large condition number of Finite-Cell operators can be linked either to a large discrepancy in the volumic fraction of the approximation elements (see Figure 20b), or to the quasilinear dependency of the polynomial basis (which tends to become the primary cause of conditioning issues for slender structures as in Figure 20a).

Preconditioning strategies, such as diagonal scaling can efficiently cure the first cause, but not the second one. De Prenter et al. proposed a Symmetric Incomplete Permuted Inverse Cholesky preconditioner (SIPIC) [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF] that can efficiently precondition for both volumic fraction and quasi-linear dependency issues. Although this preconditioning strategy is not used in this next examples of the paper, we are currently investigating its use. An initial characterization of the performances of the preconditioner is presented in Figure 21 condition number for our problem. The preconditioner is very efficient and leads to an almost constant condition number independently of the order of the approximation (for Legendre basis). The preconditioner is finally compared against classical preconditioners for two different slenderness ratios in Figure 22 (Jacobi and Incomplete Cholesky with preserved sparsity, see [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF] for the setup of the preconditioners). It can be seen that SIPIC preconditioner is to be very robust, although is seems sometimes less efficient than the incomplete Cholesky preconditioner. Nevertheless, in every case the sparsity of the SIPIC preconditioner is far better than for the incomplete Cholesky preconditioner, see Figure 23.

3D thin structures

Scordelis-Lo roof

Consider the Scordelis-Lo roof [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF] with dead load benchmark which is depicted in Figure [START_REF] Yang | An efficient integration technique for the voxel-based finite cell method[END_REF]. The geometrical and material properties are taken from [START_REF] Rank | High order finite elements for shells[END_REF]. This problem is a well known benchmark problem to assess the performances of shell and solid-shell elements. We want to treat this problem by means of the X-FEM, in order to study whether the anisotropic meshing algorithm is able to capture accurately the geometry of the roof and also its sensitivity to locking. Due to symmetries, only one fourth of the roof is modelled, and the end diaphragm is considered as rigid. The regularization parameter of the anisotropic algorithm is set to 5.10 -2 which is 40 times smaller than the thickness of the roof. The resulting mesh, which contains 18 150 nodes is depicted in Figure 25.

The computational mesh and the iso-zero of the level-set representing the geometry are shown in Figure 26: notice how coarse the mesh with respect to the geometry is.

A p convergence is conducted (p = 1 • • • 5), and both strain energy and displacement at point A (see Figure 24) are monitored. Two computational meshes are considered: one with 2 × 2 × 2 elements (as depicted in Figure 26(a)) and one with 8 × 8 × 8 elements. The vertical displacement at point A is compared to the reference value 0.00885 [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], and the strain energy to the reference value given in [START_REF] Rank | High order finite elements for shells[END_REF]: U ref = 0.003933076912. The corresponding displacement field is depicted in Figure 26

Twisted beam

Consider the Twisted Beam problem proposed in [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF], and illustrated in Figure 28. The beam is clamped on one end, and subjected to an out-of-plane loading on the opposite face (uniform tractions). This problem is known to be a severe test for assessing shell elements performance. Following [START_REF] Bouclier | An isogeometric locking-free NURBSbased solid-shell element for geometrically nonlinear analysis[END_REF], two beam thickness are considered h = 0.32 and h = 0.0032, for the so-called thick and thin cases respectively. Reference displacements at the center of the free end of the beam are given in [START_REF] Macneal | A proposed standard set of problems to test finite element accuracy[END_REF]. The problem is solved on the very coarse meshes depicted in Figure 29 by increasing the polynomial order. The convergence of the displacement of center of the free end section is presented in Figure 35, and compared to the reference solution. For the thick case, a very fast convergence is observed (for p ≥ 3), whereas the convergence is slower for the thin beam (for p ≥ 4). Yet, convergence is still very fast. For the thin beam, the use of one unique element was also tested with a very good efficiency.

Three dimensional stiffened panel

As a more realistic example, consider the stiffened panel presented in Figure 36. The general dimensions of the panel are 600 × 1000mm, and its thickness is 2.944mm. The panel is blocked along its y and z outer boundaries, and subjected to a constant pressure P on its top face. The geometry is captured by means of an anisotropic mesh composed of 64000 nodes presented in Figure 37. The regularization parameter ε was set to 10 -1 , and the bounding box of the mesh is 600 × 600 × 1000mm. Next, a p convergence is conducted, and both displacement and strain energy are compared to an overkill quadratic solid finite element computation involving 500 000 dofs.

The mechanical problem is solved on a regular grid composed of 16 elements along each side of the cube. The associated computational mesh is represented in Figure 38. Displacement and stress fields are compared for both computations (with p = 4 approximation) in Figure 39, and it can be seen that they both agree very well although the fictitious domain calculation only involves 21 000 dofs.

The corresponding convergences are depicted in Figure 40. It can be seen than both displacement and energy converge to the reference solution with no meshing burden. However, Figure 40 (b) also highlights a slower convergence rate than in the previous examples. This could be related to the geometrical accuracy that may deviate the results from the reference solution or to conditioning issues that may hinder the accuracy.

Ribbed Corner

The aim of this last example is to highlight the performance of the method for geometries involving both massive and thin areas. A ribbed corner, as depicted in Figure 41 It can be seen that the proposed approach is able to treat seamlessly both thin and thick areas of the part. The evolution of the strain energy during p refinement is presented in Figure 43. Comparison is made with a linear conforming finite element approximation with increasing number of elements, and an overkill conforming finite element solution (1 765 101 dofs) which serves as reference. It can be seen that two orders of magnitudes can be saved for a given accuracy, with respect to the direct finite element calculations. Note here that the combination of the proposed method with octree refinement near re-entrant corners would certainly improve the convergence.

Conclusion

An efficient integration method has been proposed in the context of fictitious domain methods. It relies on an edge-based error estimator and an associated metric field that drives an anisotropic mesh adaptation procedure. In 2D, it has been shown that for a given accuracy, this approach could lead to one order of magnitude less integration points compared to classical sub-grid based integration strategies. This occurs especially when the interfaces exhibit both low and large curvatures, which is common for industrial structures. For more homogeneous curvatures, the approach was demonstrated to be as least as accurate as subgrid-based integration schemes. Finally, it is important to note that this approach is top to bottom, which make it possible to discover very small features that could not be represented by the subgrid scheme if it falls below its spatial minimal resolution. It has been shown that thanks to this flexible tool, high-order fictitious domain methods could be used as solid-shell methods, handling naturally thin, solid regions and their connection. Numerical examples have demonstrated that lockingfree and computationally efficient solutions could be obtained for classical benchmarks. Still, additional work should be conducted to further improve the approach: mesh refinement strategies at the approximation level should be considered in order to capture local stress concentrations, while improving the conditioning of the operators is of primary importance if iterative solvers are used. The strategy proposed here extends naturally with tree-based approximation grids to handle the first issue, whereas preliminary results indicate that the preconditioner proposed in [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF] seems to be a promising approach.

Figure 2 :

 2 Figure 2: Interpolation of the function sin(xπ/2) cos(xπ) on two meshes composed of 100 nodes. Left: isotropic mesh, right: anisotropic mesh.

Figure 3 :

 3 Figure 3: Illustration of the anisotropic procedure for various interfaces: straight (509, 513, 517 nodes), circular (765, 787, 792 nodes) and corner (751, 759, 763 nodes)

Figure 4 :

 4 Figure 4: (a) Approximation mesh (thick lines) and integration mesh (thin lines); (b) Sliced integration mesh; (c) Resulting integration cells.

Figure 5 (Figure 5 :

 55 Figure 5: Geometrical convergence. (a): Domain ([-1, 1] 2 ) and circular interfaces (R = 0.05, 0.5, 5); (b): Convergences: octree: thick lines, anisotropic: dashed lines.

Figure 6 :

 6 Figure 6: Geometrical convergence. (a): Domain ([-1, 1] 2 ) and interfaces (R = 0.5, 5, 50); (b): Convergences: octree: thick lines, anisotropic: dashed lines.

Figure 7 :

 7 Figure 7: Infinite plate with a hole submitted to uniaxial tension: geometry and approximation mesh.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 891011 Figure 8: Typical integration meshes: octree (left) and anisotropic, ε = 10 -6 (right)

Figure 12 :

 12 Figure 12: Kinked interface: left, Mesh and domain; right, exact solution.

Figure 14 :

 14 Figure 14: Error in the energy norm. (a) with respect to dofs, (b) with respect to the number of integration points.

Figure 15 :

 15 Figure 15: Kinked interface: anisotropic mesh (639 nodes, ε = 10 -4 ).

Figure 16 :Figure 17 :

 1617 Figure 16: Quarter circle beam: problem description

Figure 18 :

 18 Figure 18: Elements containing matter for slenderness ratios of 10, 100 and 1000.

Figure 20 :

 20 Figure 20: Problem for the conditioning study. High slenderness ratio (top) and low slenderness ratio (bottom).

Figure 21 :

 21 Figure 21: Evolution of the condition number for p-refinement (p = 1 • • • 6). (a) no preconditioning; (b) SIPIC preconditioner

Figure 22 :Figure 23 :

 2223 Figure 22: Comparison of the performance of various preconditioning strategies during p-refinement (p = 1 • • • 6). (a) Slenderness = 100; (b) Slenderness = 1000

  (b) for various polynomial orders. The evolution of the displacement at point A is presented in Figure 27 (a): Even the coarsest mesh converges to the reference solution for low polynomial orders (cubic) and with less than 2000 dofs. Concerning the convergence of the relative error in the energy norm (Figure 27 (b)), very low error levels can be attained with nearly exponential

Figure 24 :

 24 Figure 24: Scordelis-Lo roof with dead load benchmark problem

Figure 25 :Figure 26 :Figure 27 :

 252627 Figure 25: Anisotropic geometrical mesh.

Figure 28 :

 28 Figure 28: Twisted beam: geometry and boundary conditions.

Figure 29 :

 29 Figure 29: Twisted beam: Geometry and approximation meshes.

Figure 30 :

 30 Figure 30: Twisted beam: Initial geometrical mesh before adaptation.

Figure 31 :

 31 Figure 31: Twisted beam: anisotropic meshes for the thick beam (ε = 10 -3 ).

Figure 32 :Figure 34 :

 3234 Figure 32: Twisted beam: anisotropic meshes for the thin beam (ε = 10 -3 ).

Figure 35 :Figure 36 :

 3536 Figure 35: Twisted beam: Convergence of the displacement field at the center of the end section. (a) Thick beam; (b) Thin beam.

Figure 37 :

 37 Figure 37: Stiffened panel: geometrical mesh (64000 nodes) and corresponding level-set iso zero.

Figure 38 :

 38 Figure 38: Stiffened panel: computational mesh.

17

 17

Figure 39 :Figure 40 :

 3940 Figure 39: Stiffened panel: comparison between the p = 4 computation (top) and the reference solution (bottom). Left: displacement, right: von Mises norm of the stress field.Note that the discrepancy between these two solutions is magnified here as the stress field is assumed to be linear on each element in the fictitious domains case (it is cubic in reality).

Figure 41 :

 41 Figure 41: Ribbed corner: (a) Geometry and boundary conditions; (b) Adapted geometrical mesh (390 473 nodes)

Figure 42 :

 42 Figure 42: Ribbed corner: Evolution of the mechanical fields during p convergence. From left to right: p = 1,2, 3, 4, 5 and 6.)

Figure 43 :

 43 Figure 43: Ribbed corner: Evolution of the strain energy during p refinement, and comparison with conforming F.E. computations.)

Acknowledgements:

The support of the ERC Grant XLS N˚291102 is gratefully acknowledged. The authors would like to thank F. de Prenter for letting an implementation of the SIPIC preconditioner available.


Xavier Roussel

Samia Boukhenouna

Sophie Rahuel-Clermont

Guy Branlant
email: guy.branlant@maem.uhp-nancy.fr

THE RATE-LIMITING STEP OF SULFIREDOXIN IS ASSOCIATED WITH THE TRANSFER OF THE γ-PHOSPHATE OF ATP TO THE SULFINIC ACID OF OVEROXIDIZED TYPICAL 2-CYS PEROXIREDOXINS

Keywords: Sulfiredoxin, rate-limiting step, Peroxiredoxin, ATP, sulfinic acid, phosphoryl anhydride intermediate Research

The rate-limiting step of sulfiredoxin is associated with the transfer of the γ-phosphate of ATP to the sulfinic acid of overoxidized typical 2-Cys peroxiredoxins

Xavier Roussel, Samia Boukhenouna, Sophie Rahuel-Clermont, Guy Branlant

Introduction

Peroxiredoxins belong to the Cys peroxidase family that reduces hydroperoxides including hydrogen peroxide (H 2 O 2 ) and peroxynitrite. Members of the typical eukaryotic two-Cys subfamily contain two essential cysteines. The first, termed the peroxidatic cysteine (C P , i.e. Cys48 in the case of the Saccharomyces cerevisiae 2-Cys-Prx Tsa1 and referred to as Prx in the following text), reduces hydrogen peroxide and alkyl hydroperoxide to water and alcohol, and is oxidized to a sulfenic acid (SOH). The second, the resolving cysteine (C R , i.e. Cys171) of the other Prx subunit within the homodimer attacks the sulfenic acid intermediate to form a disulfide [START_REF] Kang | Characterization of a mammalian peroxiredoxin that contains one conserved cysteine[END_REF][START_REF] Kong | A novel peroxiredoxin of the plant Sedum lineare is a homologue of Escherichia coli bacterioferritin co-migratory protein (Bcp)[END_REF][START_REF] Baker | Catalytic mechanism of thiol peroxidase from Escherichia coli. Sulfenic acid formation and overoxidation of essential CYS61[END_REF][START_REF] Wood | Structure, mechanism and regulation of peroxiredoxins[END_REF]. This disulfide bond is reduced by thioredoxin (Trx). The peroxidatic cysteine is sensitive to overoxidation by H 2 O 2 to sulfinic acid (SOOH) [START_REF] Wood | Structure, mechanism and regulation of peroxiredoxins[END_REF]. This overoxidation, which inactivates the enzyme, is reversed by the enzyme Sulfiredoxin (Srx) which reduces the sulfinic acid back to the sulfenic acid state [START_REF] Biteau | ATP-dependent reduction of cysteinesulphinic acid by S. cerevisiae sulphiredoxin[END_REF]. Such a reversible peroxide-sensitive switch has been demonstrated to regulate the action of certain transcription factors targeting antioxidant genes in Schizosaccharomyces pombe [START_REF] Vivancos | A cysteine-sulfinic acid in peroxiredoxin regulates H 2 O 2 -sensing by the antioxidant Pap1 pathway[END_REF].

Recently, it was shown that reduction of PrxSO 2 by Srx involved a covalent intermediate linked by a thiolsulfinate bond between the essential cysteine of Srx (referred to as the catalytic cysteine and numbered 84 in the S. cerevisiae Srx and 99 in the human enzyme) and the C P of Prx [START_REF] Jonsson | Identification of intact protein thiosulfinate intermediate in the reduction of cysteine sulfinic acid in peroxiredoxin by human sulfiredoxin[END_REF][START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. Formation of the thiolsulfinate species requires the activation of the sulfinic acid group of PrxSO 2 by transfer of the γphosphate of ATP in the presence of Mg 2+ , leading to a phosphoryl sulfinyl anhydride intermediate (Fig. 1). Reduction of the PrxSO-SSrx thiolsulfinate intermediate is completed by Trx via the formation of an intra-disulfide oxidized Srx intermediate in the case of S. cerevisiae Srx [START_REF] Roussel | Catalytic mechanism of Sulfiredoxin from Saccharomyces cerevisiae passes through an oxidized disulfide sulfiredoxin intermediate that is reduced by thioredoxin[END_REF].

Intriguingly, the catalytic rate constant of the Srx-catalyzed reaction was found to be low with a turnover of less than 2 min -1 [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF][START_REF] Chang | Characterization of Mammalian Sulfiredoxin and Its Reactivation of Hyperoxidized Peroxiredoxin through Reduction of Cysteine Sulfinic Acid in the Active Site to Cysteine[END_REF], that is limited by a step occurring before or associated with the formation of the thiolsulfinate [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. This low rate is consistent with a role of Srx in regulation of H 2 O 2dependent signaling processes since it would enable hydrogen peroxide to convey its cellular message [START_REF] Hall | Typical 2-Cys peroxiredoxins--structures, mechanisms and functions[END_REF]. This raises the question of the nature of the rate-limiting process that leads to the thiolsulfinate species. Based on rapid quench/mass spectrometry/ 18 O exchange coupled technique, a kinetic scheme was put forward for the human Srx, in which a rapid formation of the sulfinyl phosphoryl intermediate is followed by a rate-limiting attack of Srx catalytic Cys (rate constant 1.4 min -1 ) to give the thiolsulfinate intermediate, in competition with hydrolysis (rate constant 0.35 min -1 [START_REF] Jonsson | Identification of intact protein thiosulfinate intermediate in the reduction of cysteine sulfinic acid in peroxiredoxin by human sulfiredoxin[END_REF]). Thus, if the proposed scheme is right, the anhydride intermediate would accumulate.

To get insight into the molecular/ mechanistic features responsible for these unexpected slow kinetics, we used both real-time enzymatic coupled systems to determine the rates of ADP and phosphate release for wild-type, C84S and C84A Srxs from S. cerevisiae, and pH-profile analysis of the rate of the *Manuscript text (starting from the introduction) Click here to view linked References Srx-catalyzed reaction. Taken together, our data support: i) an overall rate-limiting step for the wild-type enzyme associated with the chemical process of the ATP-dependent activation of the Prx sulfinic acid that implies a non accumulation of the sulfinyl phosphoryl anhydride intermediate. Thus, as soon as the anhydride forms, either it reacts efficiently with the catalytic Cys84 of Srx or hydrolyzes when the Cys84 is substituted by an Ala residue, and; ii) a rate of formation of the phosphorylated Prx intermediate that depends on two pK apps of 6.2 and 7.5. The first pK app corresponds to the γ-phosphate of ATP and the second to an acid-base catalyst which would be involved in the ADP release.

Material and methods

Material

Tris was obtained from VWR, KCl and MgCl 2 from Merck, NADH and ATP-γ-S from Roche, and ATP, 7methylguanosine ( m7 Guo), phosphoenolpyruvate (PEP), lactate dehydrogenase from rabbit muscle (LDH) and pyruvate kinase from rabbit muscle (PK) from Sigma-Aldrich.

Recombinant purine nucleoside phosphorylase from Escherichia coli (PNP), N-terminally His 6tagged fusion of S. cerevisiae 2-Cys-Prx Tsa1 (referred to as Prx) and S. cerevisiae Srx were prepared following the experimental procedures described previously [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. Oxidation of Cys C P into the sulfinic acid state for both wild-type and C171A Prxs to form PrxSO 2 , and determination of the steady-state rate of the Srx-catalyzed reaction using the E. coli Trx/NADPH thioredoxin reductase (NTR) system with NADPH, were performed as described previously [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF].

LDH/PK/NADH-coupled assay

The kinetics of ADP release were followed under single turn-over conditions by the decrease of the absorbance at 340 nm due to the oxidation of NADH by the PK/LDH/NADH system in the presence of PEP, measured at 30 °C on an SX18MV-R stopped-flow apparatus (Applied PhotoPhysics) fitted for absorbance measurements. One syringe contained 1 mM ATP, 3 mM MgCl 2 , 1 mM PEP, and 400 µM NADH in buffer TK (50 mM TrisHCl, 100 mM KCl, pH 7) and the other contained 80 µM C171A PrxSO 2 , 20 U/ml PK, 60 U/ml LDH, and 10 µM Srx (final concentrations after mixing). For each condition, a blank time course was recorded in the absence of Srx and the data were deducted from the corresponding data obtained in the presence of Srx. An average of three runs was recorded for each set of conditions. For wild-type Srx with C171A PrxSO 2 , the rate constant k obs was obtained by fitting blank-corrected absorbance traces against equation [START_REF] Kang | Characterization of a mammalian peroxiredoxin that contains one conserved cysteine[END_REF], in which c represents the end point and a, the amplitude of the signal. For other conditions, the rate constant k obs was deduced by linear regression.

(1)

The kinetic assay was validated under steady-state conditions using wild-type and C171A PrxSO 2 (which behaves similarly) in the presence of Trx as reducing agent [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF].

PNP-coupled assay

The kinetics of Pi release were followed under single turn-over conditions as described previously [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. Briefly, the decrease of the emission fluorescence intensity associated with the phosphorolysis of m7 Guanosine ( m7 Guo) catalyzed by PNP was recorded at 30 °C on an SX18MV-R stopped-flow apparatus fitted for fluorescence measurements, with the excitation wavelength set at 305 nm, and the emitted light collected above 455 nm using a cutoff filter. One syringe contained 1 mM ATP, 1 mM MgCl 2 , 300 µM m7 Guo, and 5 µM PNP in buffer TK, and the other contained 80 µM C171A PrxSO 2 , 20 µM PNP, and 10 µM Srx (final concentrations after mixing). For each condition, a blank time course was recorded in the absence of Srx and the data were deducted from the corresponding data obtained in the presence of Srx.

For wild-type Srx with C171A PrxSO 2 , the rate constant k obs was obtained by fitting fluorescence traces against equation [START_REF] Kang | Characterization of a mammalian peroxiredoxin that contains one conserved cysteine[END_REF]. For other conditions, the rate constant k obs was obtained by linear regression. For C84A Srx, the rate constant k obs corresponded to the initial rate constant. The fluorescence signal was calibrated independently against P i concentration.

pH dependence of the rate-limiting step of the Srx-catalyzed reaction

The pH-profile of the rate-limiting step of the reaction was determined by following the decrease of intrinsic tryptophan fluorescence intensity of the wild-type Prx substrate upon transitioning from the SO 2 form to the disulfide form. As shown previously for the S. cerevisiae wild-type enzyme in the absence of added reductant [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF], the reaction proceeds until the oxidized disulfide form of Prx is formed, whose intrinsic fluorescence signal is quenched by approximately 80% compared to the PrxSO 2 substrate. The decrease of the emission fluorescence intensity was recorded at 30°C on an SX18MV-R stopped-flow apparatus fitted for fluorescence measurements, with excitation wavelength set at 295 nm, and emitted light collected above 320 nm using a cutoff filter. One syringe contained 1 mM ATP, 1 mM MgCl 2 , and 40 µM PrxSO 2 , and the other contained 5 µM Srx (final concentrations after mixing). Each measurement

c ae y t k obs    -4 -
        
was carried out in 25 mM Tricine, 100 mM imidazole, and 75 mM acetate buffer at an constant ionic strength of 0.1 M over a pH range of 5-9 [START_REF] Ellis | Buffers of constant ionic strength for studying pH-dependent processes[END_REF]. The ionic strength was verified by conductimetry and adjusted to 0.1 M when necessary using KCl, after calibration against KCl solutions.

For each condition, a blank time course was recorded in the absence of Srx and the data were deducted from the corresponding data obtained in the presence of Srx. The rate constant k obs was obtained by fitting the fluorescence traces against equation [START_REF] Kang | Characterization of a mammalian peroxiredoxin that contains one conserved cysteine[END_REF]. The pH profile of the deduced rate constants was then fitted to equation ( 2), in which pK 1 and pK 2 represent the measured apparent pK values, k obsmax the maximum rate constant and γk obsmax the rate constant at high pH.

(2)

Results

Characterization of the rate-limiting step

To identify which step(s) in the reaction catalyzed by the Srx from S. cerevisiae contributes to the rate of formation of the thiolsulfinate intermediate, the LDH/PK/NADH coupled assay was used to monitor the quantity of ADP released in the first step. In this assay, phosphorylation of ADP by PK in the presence of PEP leads to production of pyruvate, which is then reduced by LDH with the consumption of NADH. Thus, the rate of ADP formation which can be monitored indirectly via the decrease in absorbance of NADH was determined under single turn-over conditionsthat is, in the presence of an 8fold excess of C171A PrxSO 2 relative to wild-type Srx, and in the absence of reductant.

A control reaction lacking wild-type Srx produced a slow linear signal with rate lower than 0.05 min -1 which was subtracted from traces collected under the same conditions in the presence of wild-type Srx.

The best fit to the resulting progress curves was obtained with a monoexponential process, yielding a first order rate constant of 2.8 ± 0.5 min -1 which was independent of wild-type Srx concentration (Fig. 2a). This rate constant is close to the rate of Pi release determined under single turnover conditions (1.8 ± 0.2 min -1

) and to the steady-state rate constant of the reaction measured either by the NADPH/Trx/NTR coupled assay (1.9 ± 0.2 min -1 ) [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF] or by following the rate of Pi release by the PNP-coupled assay (1.7 ± 0.2 min -1 ) in the presence of excess of Trx [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. The amplitude of the exponential process was proportional to the concentration of wild-type Srx.

-5 -

Rates of release of ADP and phosphate by C84A and C84S Srxs

The rates of release of ADP and phosphate by C84A and C84S Srxs were determined during reaction with C171A PrxSO 2 in the absence of reductant. As described previously, the release of ADP was monitored by the LDH/PK/NADH coupled assay while that of phosphate was followed by the PNP coupled assay. For C84A Srx, the time course of ADP release was linear, yielding a rate constant of 0.6 min -1 (Fig. 2a). We observed a decrease in the rate of phosphate release with time, which resulted in a non-linear time course. Such a behavior remains to be explained. The rate constant for phosphate production in this case was deduced from the initial slope of the curve (Fig. 2b), giving a value of 0.7 min -1 , which is similar to that of ADP release. The same experiment with C84S Srx yielded an essentially linear time course of ADP release with a rate constant of 0.1 min -1 (Fig. 2a), but no phosphate production was seen (Fig. 2b). A control experiment was carried out under the same conditions using reduced Prx instead of the overoxidized substrate and C84A Srx to evaluate the rate of background ADP and phosphate production. A rate constant of 0.05 min -1 was measured for ADP release while release of free phosphate was not observed. Thus, the phosphate which is released in the presence of C84A Srx and the substrate PrxSO 2 can only derive from hydrolysis of the sulfinyl phosphoryl anhydride species formed on the overoxidized substrate. Indeed, an alternative origin for the phosphate via an ATPase activity of C84A Srx is unlikely, given that the background rate of release of ADP is only 0.05 min -1 , and background production of phosphate is undetectable for this mutated Srx.

pH dependence of the rate-limiting step of the Srx-catalyzed reaction

We did not determine the pH-rate profiles using the PNP-or LDH/PK/NADH-coupled assays due to a potential dependence of the assays themselves on pH. Instead, we used as a kinetic probe the attenuation of intrinsic tryptophan fluorescence of wild-type Prx upon formation of the oxidized disulfide species PrxC P S-SC R Prx, which is rapid and follows the rate of release of the reaction product PrxSOH [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF], whose rate is itself limited by that of the activation of the Prx sulfinic acid. The experiments were carried out under sub-saturating concentrations of PrxSO 2 (2 K M ) to obtain a sufficient signal to noise ratio, the signal being proportional to the Srx concentration, relative to the intrinsic fluorescence of the PrxSO 2 . As shown in Fig. 3a, the rate of reaction exhibits a bell-shaped dependence on pH, which indicates that at least two groups are required to be in the proper ionization state for maximal catalytic activity at pH 7. One of these groups, characterized by a pK app of 6.2 ± 0.2, must be deprotonated, whereas the other with a pK app of 7.5 ± 0.2, must be protonated. When γ-S-ATP, whose K M value is similar to that of ATP (Table 1), was used instead of ATP, the first pK app value shifted from 6.2 to 5.5, while the second was determined to be of 8.0 ± 0.2 (Fig. 3b). Moreover, a two-fold decrease in the k obs max was observed (Table 1). Using Mn 2+ instead of Mg 2+ led i) the first pK app value to shift from 6.2 to 5.5, and the second from 7.5 to 7.8; and ii) to similar k obs max (Fig. 4c) and K M values for Mn 2+ and Mg 2+ (Table 1). In principle, these two pK apps can correspond to ionizable groups on amino acids of Prx or Srx within the quaternary Srx/Prx/ATP/Mg 2+ complex or to the γ-phosphate of ATP. Because the concentration of PrxSO 2 is not fully saturating, the k obs vs. pH curve could also include a PrxSO 2 k obs /K M contribution.

Discussion

In the present study, the rates of release of phosphate and ADP by wild-type Srxs from S. cerevisiae have been determined under single turnover conditions. In the case of wild-type Srx with C171A PrxSO 2 , the rate of ADP release is similar to that of phosphate release, both determined under single turnover conditions (2.8 min -1 vs. 1.8 min -1 ) [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF]. As a thiolsulfinate intermediate between Prx and wild-type Srx was shown previously to be formed at a rate of 1.7 min -1 under steady-state conditions [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF], this finding supports an overall rate-limiting step for the wild type that includes all steps up to the release of ADP and thus before the formation of the thiolsulfinate intermediate.

In a previous study, a stoichiometry of one mole of phosphate released per mole of wild-type Srx was observed upon formation of the thiolsulfinate intermediate in the absence of added reductant [START_REF] Roussel | Evidence for the formation of a covalent thiosulfinate intermediate with peroxiredoxin in the catalytic mechanism of sulfiredoxin[END_REF].

Combined with the above result, this demonstrates that hydrolysis of the sulfinyl phosphate does not occur, nor is detectable under the experimental conditions used. This finding implies that each sulfinyl phosphate anhydride Prx intermediate which is formed proceeds to the thiolsulfinate Prx/Srx species at a much higher intrinsic rate than the 1.4 min -1 reported [START_REF] Jonsson | Identification of intact protein thiosulfinate intermediate in the reduction of cysteine sulfinic acid in peroxiredoxin by human sulfiredoxin[END_REF] and does not accumulate. The reason why the Prx intermediate does not hydrolyze but forms a thiolsulfinate Prx/Srx intermediate is the higher chemical reactivity of the Srx Cys84 (whose pK app is 7.4 1 -as also shown for the human Srx [START_REF] Chang | Characterization of Mammalian Sulfiredoxin and Its Reactivation of Hyperoxidized Peroxiredoxin through Reduction of Cysteine Sulfinic Acid in the Active Site to Cysteine[END_REF]) towards the sulfinyl phosphate anhydride intermediate, relative to a non activated water molecule, a reactivity which is further favored by an entropic effect due to the relative positioning of the anhydride intermediate and Cys84 within the Prx/Srx complex.

To get deeper into the mechanism of phosphate transfer from ATP to the sulfinate group of PrxSO 2 , we performed a kinetic analysis of ADP release rates using C84A and C84S mutants of Srx. As deduced from Fig. 2, the rates of release of ADP and phosphate with C84A Srx are in the same range and significantly higher (0.6-0.7 min -1 ) than the background. Therefore, formation of a phosphoryl-Cys intermediate on the Cys84 of wild-type Srx can be excluded, in contrast to that suggested by Jeong and co-workers [START_REF] Jeong | Molecular mechanism of the reduction of cysteine sulfinic acid of peroxiredoxin to cysteine by mammalian sulfiredoxin[END_REF]. It is also in accord with structural data obtained on the human enzyme [START_REF] Jonsson | Reduction of cysteine sulfinic acid in peroxiredoxin by sulfiredoxin proceeds directly through a sulfinic phosphoryl ester intermediate[END_REF]. The fact that both rates are similar and the time course, in particular of ADP release, is linear indicates no accumulation of the sulfinyl phosphoryl anhydride intermediate. Therefore, the 'intrinsic' rate of hydrolysis of this intermediate is likely significantly higher than 0.6 min -1 . Given the difference of pK a between sulfinic and acetic acid (1.9 vs. 4.3, respectively), this result is consistent with the rate of hydrolysis of a phosphoryl carboxylic anhydride like acetyl phosphate, (0.004 and 0.01 min -1 at pH 7.7 in the absence and presence of 100 mM of Mg 2+ , respectively [START_REF] Koshland | Effect of Catalysts on the Hydrolysis of Acetyl Phosphate. Nucleophilic Displacement Mechanisms in Enzymatic Reactions1[END_REF]). Due to the regeneration of the substrate PrxSO 2 , a futile cycle occurs that consumes ATP with linear, steady state kinetics, and an amplitude that is only limited by ATP and Mg 2+ concentration, contrary to the wild-type mechanism that corresponds to single-turnover kinetics, with amplitude limited by the concentration of Srx.

The decrease in the rates of ADP and phosphate release for C84A Srx relative to wild type when using C171A PrxSO 2 as a substrate could be caused by a change in the environment near the γ-phosphate of ATP due to the substitution of Cys by Ala. In the case of C84S Srx, the rate of release of ADP is 0.1 min -1 , which is only two-fold higher than the background rate, and 28 and 6-fold lower than the rate observed with wild-type and C84A Srxs, respectively. These data suggested that for as yet undetermined structural reasons which is likely related to the high flexibility of the active site of Srx, the quaternary PrxSO 2 /ATP/Mg 2+ /C84S Srx complex is not competent to transfer the γ-phosphate of ATP to PrxSO 2 .

The fact that no phosphate release was observed supports this interpretation.

The pH dependence of the rate of formation of the sulfinyl phosphoryl anhydride shows a bellshaped profile from which two pK apps of 6.2 and 7.5 were extracted. Based on the above demonstration that phosphate transfer does not depend on Cys84, the first ionization of pK app 6.2 may be attributed to either the sulfinate group of PrxSO 2 , or to the γ-phosphate group of ATP. Attribution of the ionization to the γ-phosphate was established from the observation that the first pK app decreased when ATP-γ-S and Mn 2+ were used instead of ATP [START_REF] Pecoraro | Stability constants of Mg 2+ and Cd 2+ complexes of adenine nucleotides and thionucleotides and rate constants for formation and dissociation of MgATP and MgADP[END_REF][START_REF] Frey | Bond order and charge localization in nucleoside phosphorothioates[END_REF] and Mg 2+ [START_REF] Sigel | Isomeric equilibria in complexes of adenosine 5'-triphosphate with divalent metal ions. Solution structures of M(ATP) 2 -complexes[END_REF], respectively. Taken together, the data clearly demonstrate: i) a rate-limiting step associated with the chemical process leading to the transfer of the γphosphate from ATP to the sulfinic acid; and ii) the γ-phosphate of ATP as the group with a pK app of 6.2, and thus excludes the sulfinic acid whose pK app in solution is ~2 [START_REF] Reddie | Expanding the functional diversity of proteins through cysteine oxidation[END_REF].

The nature of the second pK app remains to be determined. The difficulty with this measurement stems partly from the overall low rate of turnover of Srx, which precludes kinetic comparison of the effects of different mutations introduced into the active site. The transfer of the γ-phosphate involves the release of ADP which is likely to be assisted by an acid catalyst. A candidate residue deduced from the inspection of the crystal structure of human Srx in complex with ATP-Mg 2+ is His100 [START_REF] Jonsson | Reduction of cysteine sulfinic acid in peroxiredoxin by sulfiredoxin proceeds directly through a sulfinic phosphoryl ester intermediate[END_REF], which corresponds to His85 in the S. cerevisiae Srx. Indeed, substitution of this His by Asn either in the S. cerevisiae 1 or human [START_REF] Jeong | Molecular mechanism of the reduction of cysteine sulfinic acid of peroxiredoxin to cysteine by mammalian sulfiredoxin[END_REF] enzymes leads to a complete loss of activity, with a k cat value approximately equal to that of background (e.g. a 20-fold decrease at least). -1 -

Figure 1 .

 1 Figure 1. Catalytic mechanism of the formation of the thiolsulfinate intermediate of Srx. P: Prx; S: Srx.

Figure 2 .

 2 Figure 2. Single turnover kinetics of the reaction catalyzed by Srx with C171A PrxSO 2 . Reactions of 80 µM C171A PrxSO 2 with 10 µM wild-type (black), C84S (grey) and C84A (light grey) Srxs in the presence of 1 mM ATP/MgCl 2 were carried out in buffer TK at 30 °C. (a) The reaction was monitored by ADP release by the LDH/PK/NADH coupled assay (see section 2.2). (b) The reaction was monitored by the PNP coupled assay (see section 2.3). The progress curves are corrected from data collected in the absence of Srx.

Figure 3 .

 3 Figure 3. Representative pH dependence of the rate-limiting step of the Srx-catalyzed reaction. Reaction of 40 µM wild-type PrxSO 2 and 5 µM wild-type Srx was followed at 30°C by the attenuation of the intrinsic tryptophan fluorescence of PrxSO 2 (see section 2.4) in the presence of 1 mM ATP, 1 mM MgCl 2 (a), 1 mM ATP-γ-S and 1 mM MgCl 2 (b) or 1 mM ATP, 1 mM MnCl 2 (c). Each measurement was carried out in 25 mM tricine, 100 mM imidazole, and 75 mM acetate buffer at a constant ionic strength of 0.1 M. Apparent pK 1 and pK 2 values of 6.2 ± 0.2 and 7.5 ± 0.2 (a); 5.5 ± 0.1 and 8.0 ± 0.2 (b); and 5.5 ± 0.1 and 7.8 ± 0.2 (c) were deduced from the best fit (solid line) of the experimental points (•) against equation (2). Results are expressed as the average of two independent determinations (± SD).

Table 1 . Kinetic parameters of wild-type Srx for the co substrates ATP and Mg 2+ , and for the co substrate analogs γ-thio-ATP and Mn 2+ .

 1 The Michaelis constants were determined from activity measurement performed using the Trx/NTR coupled assay in TK buffer pH 7 at 30°C by following the consumption of NADPH at 340 nm. Reaction of 50 µM PrxSO 2 , 200 µM NADPH, 50 µM Trx, 0.5 µM NTR and 1 mM MgCl 2 or 1 mM ATP with variable concentrations of the indicated co substrate, were initiated by adding 10 µM Srx.

		ATP	γ-thio-ATP	MgCl 2	MnCl 2
	K M (µM)	125 ± 50	130 ± 30	250 ± 60	230 ± 70
	k obs max (min -1 )	1.83 ± 0.18	0.96 ± 0.04	1.82 ± 0.13	1.75 ± 0.15

Table 1 Figure1 Click here to download high resolution image Figure2 Click here to download high resolution image Figure3 Click here to download high resolution image

 1

Acknowledgements

This work was supported in part by the CNRS, the University of Nancy I, the Institut Fédératif de Recherche 111 Bioingénierie, and by the French Agence Nationale de la Recherche (program ANR-06-BLAN-0369). X. Roussel and S. Boukhenouma was and is supported, respectively by the French Research Ministry. We thank I. Artaud and K. Weissman for helpful discussions and A. Kriznik and J.

Charbonnel for excellent technical assistance.

From a chemical point of view, the rate of the phosphotransferase activity was expected to be rapid. Indeed, a sulfinic acid has a pK app of ~2 and is a better nucleophile than a carboxylate group [START_REF] Reddie | Expanding the functional diversity of proteins through cysteine oxidation[END_REF].

However, the measured k obs max value is low (2 min -1 ). The structural/molecular factors which would be responsible for this low k obs max value remains to be elucidated. Nevertheless, the rate of formation of the sulfinyl phosphoryl anhydride intermediate, although low, is sufficient to repair the overoxidized Prx as shown by in vivo experiments [START_REF] Biteau | ATP-dependent reduction of cysteinesulphinic acid by S. cerevisiae sulphiredoxin[END_REF][START_REF] Chang | Characterization of Mammalian Sulfiredoxin and Its Reactivation of Hyperoxidized Peroxiredoxin through Reduction of Cysteine Sulfinic Acid in the Active Site to Cysteine[END_REF]. This fact could explain why no additional evolutionary pressure has occurred to force the sulfinic acid to reveal its 'intrinsic' chemical reactivity.


Célestin Matte

Mathieu Cunche

Spread of MAC address randomization studied using locally administered MAC addresses use historic

Keywords: Wi-Fi, tracking, privacy, random MAC address, locally administered bit

In this document, we study the spread of random MAC addresses in Wi-Fi service discovery. To do so, we gather several datasets of probe requests from 2013-2017 and leverage an indicator of such addresses: the Locally Administered bit. We observe a trend of global increase in absolute random addresses use, even though the per-frame count these addresses is still low.

Introduction

As a defense against the possibility of tracking, various vendors of Wi-Fi-enabled devices started implementing MAC address randomization during Wi-Fi service discovery [START_REF] Matte | Wi-Fi Tracking: Fingerprinting Attacks and Counter-Measures[END_REF]. In this document, we study the spread of such random addresses in the wild.

In order to communicate, devices address each other on the MAC layer using a 6-byte globallyunique identier called the MAC address. As seen in gure 1, the rst three bytes (prex) of this address are an Organizationally Unique Identier (OUI) which has to be bought by vendors from the IEEE Registration Authority in order to be used, so as to ensure the global uniqueness of MAC addresses. The last three bytes are the Network Interface Controller (NIC). One very specic bit of the MAC address is the seventh bit of the rst byte of the OUI: the Locally Administered bit (LA bit). If set to 1, it indicates that the MAC address has been changed by the administrator of the device and is not guaranteed to be unique. A prex whose LA bit is set to 1 is called a Company ID (CID), and is automatically purchased by manufacturers along with their OUI equivalent [START_REF]IEEE. Guidelines for use organizationally unique identier (oui) and company[END_REF]. It is possible for companies to buy private OUIs, which are not publicly tied to the company's name [START_REF]IEEE. Guidelines for use organizationally unique identier (oui) and company[END_REF].

It is unclear whether MAC address randomization during active scanning should set the LA bit to 1, as, to our knowledge, no document mentions this case explicitly [START_REF]IEEE. Guidelines for use organizationally unique identier (oui) and company[END_REF]. While all rmware implementing randomization do not respect this behaviour, it constitutes a good indicator of MAC address randomization use among time. Apart from randomization, this bit is only set when users manually change their MAC address, or for rarely-used protocols (e.g. Wi-Fi Direct).

We note, however, that some common devices (Nintendo devices), operate in P2P using an address with LA bit set to 1 [START_REF] Martin | A study of MAC address randomization in mobile devices and when it fails[END_REF]. Numbers presented later in this section indicate that, at a time when randomization was not integrated by vendors, LA bit use was anecdotal.

As no public study following the evolution of MAC address randomization exists, we compare LA bit use in datasets from dierent times and locations in table 1. This table lists LA bit use and the fraction of unallocated OUIs in randomized addresses in dierent datasets.

datasets

Used datasets are summarized in gure 2. We distinguish the situation of these datasets into several cases: hotspot: the recording matches what one would obtain by recording trac at a public hotspot, i.e., users stay for some time (from minutes to hours) and the turn-over is high, To give more details, we compare the fraction of vendors (according to OUIs) in the dierent datasets we have access to (in terms of number of non-random MAC addresses, identied by their LA bit) in gure 3. This gives an overview of the distribution of devices' vendors in the dierent datasets. Apple arrives rst in all datasets, as the company manufactures all iOS devices, unlike Google and its Android OS. Large dierences of appearance of a given vendor across several datasets can be attributed to their dates. Lower scores in all popular manufacturers for the Madeira can be explained by a hardly quantiable number of random addresses not setting the LA bit to 1, and using already-attributed OUIs 2 . We're unsure about this source of incoherent addresses: it could be explained by some unidentied device models performing MAC address randomization in a yet-unobserved manner, or by an attacker deliberately spamming the systems with fake addresses 3 .

RR n°9142

Methodology

We compute the amount of devices setting their LA bit to 1 while emitting probe request frames.

As we cannot reliably group probe requests from devices using MAC address randomization, we cannot obtain an exact per-device count. As a solution to this issue, we use both a per-MAC-address count and a per-probe-request count to compute the results. The former gives information on the actual number of detected MAC addresses, while the latter more accurately matches a per-device count. We could not compute the per-probe requests counts for 2 datasets because we did not have access to them. We must note that the per-MAC-address counts are expectedly highly variable: as devices using randomization use multiple MAC addresses, a small number of devices can drastically increase the amount of received randomized addresses. This is especially true if devices are monitored for a long period. For example, compare both counts for the Middleware2014 dataset in the result (see next section). 1 The Sapienza contains all these cases in dierent capture les. 2 The large number of such addresses using OUIs of small companies not manufacturing any mobile devices tends to make us believe that these are not genuine OUIs. 3 The attacker hypothesis seems more credible when we consider the fact that almost no MAC addresses in this dataset (0.05%) use both an unregistered OUI and a LA bit not set to 1. Represented manufacturers are those for which at least one dataset had 3% of its MAC addresses belonging to it. [START_REF] Matte | Wi-Fi Tracking: Fingerprinting Attacks and Counter-Measures[END_REF] Results

Inria

Results are presented in table 1.

This table indicates a clear evolution of LA bit use among time. Anecdotal in 2013, it impacted an important portion of probe requests in the beginning of 2016. Surprisingly, a very high amount of randomized addresses in most datasets use an unallocated OUI. Vendors are supposed to use registered OUIs when setting local addresses [START_REF]IEEE. Guidelines for use organizationally unique identier (oui) and company[END_REF] When considering the request count in recent datasets, results depend on the type of dataset. In a local AP conguration (Glimps2015 dataset), the fraction is close to the one obtained with a per-MAC address count, as most devices are detected for a brief amount of time.

In a hotspot conguration (Belgium dataset), devices are detected for a long amount of time, resulting in a lot of address changes. As a consequence, while a lot of dierent addresses are captured (almost 50% of which are random), it amounts for a small fraction of probe requests (less than 3%). The latter case suggests that the fraction of devices using randomization is actually low in this dataset.

Conclusion

This study reveals a trend of increasing adoption of MAC address randomization. Most of these random addresses use an unregistered OUI (more than 95% in all recent datasets). While the fraction of detected random addresses is high in recent datasets (2016), a per-frame count of the use of these addresses suggests that the number of devices using MAC address randomization is still pretty low (less than 3%).

As a future work, it would be interesting to compute this trend on a given dataset running over a long period of time, to see the trend in a xed environment.

Figure 1 :

 1 Figure 1: MAC address format.

Figure 2 :

 2 Figure 2: Used datasets

Figure 3 :

 3 Figure 3: Fraction of non-random MAC addresses belonging to most-spread manufacturers.

Table 1 :

 1 Fraction of MAC addresses having a Locally Administered bit set to 1 over the total number of MAC addresses, in dierent datasets. LA bit % columns indicate the fraction of MAC addresses having their LA bit set to 1. The Unalloc. column indicates fraction of these random addresses also using an unallocated OUI. Results are displayed by MAC addresses counts, and by probe requests count.

		Dataset		Results Per MAC addr. Per probe requests
	Time	Name	MAC addr.	Probe req. LA bit % Unalloc. LA bit % Unalloc.
	13.02-13.05 Sapienza	160 000	8 000 000	0.2%	33%	0.2%	13.5%
	14.12	Middleware2014	900	140 000	47%	99.8%	1.5%	99.8%
	15.10	Lab	1 300	120 000	14%	100%	1.7%	100%
	15.10-15.11 Train station	9 700	110 000	23%	97.2%	10.0%	89.1%
	15.12	Glimps2015	83 000	120 000	66.2%	99.0%	57.7%	98.9%
	16.01-16.02 Belgium	3 700	200 000	48.8%	99.3%	2.8%	99.3%
	15.01-16.12 Martin	2 600 000	66 000 000	53.8%	95.5%
	15.12-17.06 Madeira	13 000 000 300 000 000	99.8%	99.4%

. In fact, the equivalent CID of an OUI is automatically bought. This might indicate that most of these received probe requests come from implementations of MAC address randomization not respecting the 802.11 norm regarding MAC address use, such as manually-installed randomization applications. An example of this would be a device manually congured using the macchanger tool with its -b option (burned-in-address option, i.e., do-not-set-LA-bit option). The per-probe request counts indicate that these random addresses actually account for a very small fraction of devices, except in the Glimps2015 dataset. This exception is not really surprising, considering the fact that devices are seen for a short time in this dataset.In the Martin dataset, a progression of several Android-related CIDs used for MAC address randomization raises the fraction of registered OUI to 4.4%. However, the release of the iOS 10 OS around the same period (2016-09) also plays in favor of the former trend of unallocated OUI use for randomization, as iOS 10 devices use unregistered OUIs when using random addresses[START_REF] Martin | A study of MAC address randomization in mobile devices and when it fails[END_REF].In 2012, Musa et al. remarked a similar but slightly dierent trend: in a 9-month public street deployment, 14% of the observed unique MAC addresses (over 60 000) used an unallocated OUI, whose LA bit is not set[START_REF] Musa | Tracking unmodied smartphones using Wi-Fi monitors[END_REF]. We did not observe this trend in any of the datasets we have access to: this ratio barely reaches 1% in one of them, and lies between 0.0 and 0.1% in the other ones.[START_REF] Matte | Wi-Fi Tracking: Fingerprinting Attacks and Counter-Measures[END_REF] This document, however, does not explicitly mention systematic randomization of addresses during active scanning. To our knowledge, no ocial document does so.RR n°9142