Gabriel Caloz

Monique Dauge

VICTOR P ÉRON Erwan Faou

ON THE INFLUENCE OF THE GEOMETRY ON SKIN EFFECT IN ELECTROMAGNETISM

We consider the equations of electromagnetism set on a domain made of a dielectric and a conductor subdomain in a regime where the conductivity is large. Assuming smoothness for the dielectric-conductor interface, relying on recent works we prove that the solution of the Maxwell equations admits a multiscale asymptotic expansion with profile terms rapidly decaying inside the conductor. This skin effect is measured by introducing a skin depth function that turns out to depend on the mean curvature of the boundary of the conductor. We then confirm these asymptotic results by numerical experiments in various axisymmetric configurations. We also investigate numerically the case of a nonsmooth interface, namely a cylindrical conductor.

INTRODUCTION

Our interest lies in the influence of the geometry of a conducting body on the skin effect in electromagnetism. This effect describes the rapid decay of electromagnetic fields with depth inside a metallic conductor. The skin effect reflects the flow of current near the surface of a conductor. After the early work [START_REF] Rytov | Calcul du skin effect par la méthode des perturbations[END_REF], the mathematical analysis of the skin effect has been addressed more recently in several papers, [START_REF] Stephan | Solution procedures for interface problems in acoustics and electromagnetics[END_REF][START_REF] Maccamy | Solution procedures for three-dimensional eddy current problems[END_REF][START_REF] Maccamy | A skin effect approximation for eddy current problems[END_REF][START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF].

The present work is motivated by recent studies [START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF][START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF][START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF] in which authors analyze the behavior of the electromagnetic fields solution of the Maxwell equations through an asymptotic expansion for large conductivity. In particular, uniform estimates for the electromagnetic field at high conductivity are proved in [START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF], whereas in the note [START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF] a suitable skin depth function is introduced on the interface between a conductor and an insulator to generalize the classical scalar quantity. An asymptotic expansion at high conductivity for this function shows the influence of the geometry of the interface : the skin depth is larger for high conductivity when the mean curvature of the conducting body surface is larger -and here the sign of the curvature has a major influence, which means that the skin depth is larger in convex than in concave conductors.

In this paper, our aim is twofold [START_REF] Bernardi | Spectral methods for axisymmetric domains[END_REF] Present elements of derivation for asymptotic expansions near the conductorinsulator interface, (2) Illustrate by numerical computations the theoretical behavior deduced from asymptotic analysis in [START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF].

For our computations, we consider a special class of axisymmetric problems, which allows their reduction to two-dimensional scalar problems. This enables to measure the skin effect with reduced computational effort. We perform finite element computations for relevant benchmarks (cylinders and spheroids). The interface problem for the magnetic field is solved with the library Mélina [START_REF] Mélina | bibliothèque de calculs éléments finis[END_REF]. A postprocessing shows the accuracy of the asymptotic expansion, and exhibits the influence of the mean curvature of the conductor on the skin effect. Our computations also clearly display the effect of the edges of the cylindrical conductor: In this case, the decay is not exponential near edges while becoming exponential further away from edges.

The presentation of the paper proceeds as follows. In section 2, we introduce the framework: we present the Maxwell equations and the result of [START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF] about existence of solution at high conductivity and the formulation for the magnetic field. In section 3, we present the asymptotic expansion of the magnetic field in the smooth case, and the asymptotic behavior of the skin depth for high conductivity. In section 4, we restrict our considerations to axisymmetric domains and orthoradial axisymmetric data, and we present the configurations A (cylindrical), B, and C (spheroidal), chosen for computations. In section 5, we introduce a finite element discretization for the solution of the problem in the meridian domain, and we check the convergence of the discretized problem. In section 6, we present numerical simulations in all configurations to highlight the skin effect, and to exhibit the influence of the mean curvature of the conductor on the skin effect. In section 7, we perform post-treatments of our numerical computations in configurations A (cylindrical) and B (spheroidal) in order to investigate the nature of the decay of the field inside the conductor: In smooth configurations (B) the rate of the exponential decay is very close to the expected theoretical one, while in corner configurations (A) the exponential decay shows up in a region which is not very close to the corner. Theoretical aspects in this latter case are going to be studied in [START_REF] Dauge | Asymptotic expansion for the solution of a stiff transmission problem in electromagnetism with a singular interface[END_REF]. In Appendix A, we provide elements of proof for the multiscale expansion given in section 3. We expand the Maxwell operators in power series, and we obtain the equations satisfied by the magnetic and electric profiles. In Appendix B, we derive the profiles of the orthoradial component of the magnetic field in relation with the 3D asymptotic expansion in axisymmetric configurations.

FRAMEWORK

Let Ω be a piecewise smooth Lipschitz domain. We denote by Γ its boundary and by n the outer unit normal field on Γ. We denote by L 2 (Ω) the space of three-component fields square integrable on Ω.

We consider the Maxwell equations given by Faraday's and Ampère's laws in Ω:

(2.1) curl E -iωµ 0 H = 0 and curl H + (iωε 0 -σ)E = j in Ω .
Here, (E, H) represents the electromagnetic field, µ 0 is the magnetic permeability, ε 0 the electric permittivity, ω the angular frequency, j represents a current density and is supposed to belong to

H 0 (div, Ω) = {u ∈ L 2 (Ω) | div u ∈ L 2 (Ω), u • n = 0 on Γ} ,
and σ is the electric conductivity. We assume that the domain Ω is made of two (connected) subdomains Ω + and Ω -in which the coefficient σ take two different values (σ + = 0, σ -≡ σ). We denote by Σ the interface between the subdomains Ω + and Ω -. In the most part of our analysis 1 , we will assume that the external boundary Γ is contained in the dielectric region Ω + , see Figure 1.

Ω - Ω + Σ Γ FIGURE 1.
The domain Ω and its subdomains Ω + (dielectric) and Ω -(conductor).

To complement the Maxwell harmonic equations (2.1), we consider the perfectly insulating electric boundary conditions on the boundary Γ

(2.2) E • n = 0 and H × n = 0 on Γ .
Note that we could also consider the perfectly conducting electric boundary conditions E × n = 0 and H • n = 0 on Γ.

1 Except in one configuration in which the roles of the conductor and dielectric bodies are swapped.

2.1. Existence of solutions. Subsequently, we assume that the following condition on the limit problem in the dielectric part Ω + is valid:

Assumption 2.1. The angular frequency ω is not an eigenfrequency of the problem

(2.3)    curl E -iωµ 0 H = 0 and curl H + iωε 0 E = 0 in Ω + E × n = 0 and H • n = 0 on Σ (2.2) on Γ.
Hereafter, we denote by

(2.4) E 0,Ω + H 0,Ω + √ σ E 0,Ω -C j H(div,Ω) .
For convenience, we introduce in the Maxwell equations (2.1) the small parameter

(2.5) δ = ωε 0 /σ .
Hence, δ tends to 0 when σ → ∞. For σ σ 0 , we denote by (E (δ) , H (δ) ) the solution of the system (2.1) -(2.2).

Magnetic formulation.

By a standard procedure we deduce from the Maxwell system (2.1)-(2.2) the following variational formulation for the magnetic field H (δ) . The variational space is H 0 (curl, Ω):

(2.6) H 0 (curl, Ω) = {u ∈ L 2 (Ω) | curl u ∈ L 2 (Ω), u × n = 0 on Γ} ,
and the variational problem writes

Find H (δ) ∈ H 0 (curl, Ω) such that for all K ∈ H 0 (curl, Ω) (2.7) Ω 1 ε(δ) curl H (δ) • curl K -κ 2 H (δ) • K dx = Ω curl j ε(δ) • K dx ,
where we have set

(2.8) ε(δ) = 1 Ω + + (1 + i δ 2 ) 1 Ω -and κ := ω √ ε 0 µ 0 .
Assumption 2.3. We assume that the surfaces Σ (interface) and Γ (external boundary) are smooth.

Under Assumption 2.3, the magnetic field H (δ) is globally in H 1 (Ω).

MULTISCALE EXPANSION

Several works are devoted to asymptotic expansions at high conductivity of the electromagnetic field in a domain made of two subdomains when the interface is smooth: see [START_REF] Stephan | Solution procedures for interface problems in acoustics and electromagnetics[END_REF][START_REF] Maccamy | Solution procedures for three-dimensional eddy current problems[END_REF][START_REF] Maccamy | A skin effect approximation for eddy current problems[END_REF] for plane interface and eddy current approximation, and [START_REF] Haddar | Generalized impedance boundary conditions for scattering problems from strongly absorbing obstacles: the case of Maxwell's equations[END_REF][START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF][START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF] for a three-dimensional model of skin effect in electromagnetism.

In this section we recall the results presented in [START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF] on the behavior of the skin depth function. This relies on asymptotic expansions for the electromagnetic field analyzed in the PhD thesis [START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]. For the sake of completeness we give in Appendix A elements of proofs.

Here, we assume that j is smooth and that Assumptions 2.1 and 2.3 hold. Let U -be a tubular neighborhood of the surface Σ in the conductor part Ω -, see Figure 2. We denote by (y α , y 3 ) a local normal coordinate system to the surface Σ in U -: Here, y α , α = 1, 2, are tangential coordinates on Σ and y 3 is the normal coordinate to Σ, cf. [START_REF] Naghdi | Foundations of elastic shell theory[END_REF][START_REF] Faou | Elasticity on a thin shell: formal series solution[END_REF].

Ω - h Σ Σ h n U - FIGURE 2.
A tubular neighbourhood of the surface Σ In view of our numerical computations, we concentrate on the magnetic field H (δ) solution of (2.7). It is denoted by H + (δ) in the dielectric part Ω + , and by H - (δ) in the conducting part Ω -. Both parts exhibit series expansions in powers of δ:

H + (δ) (x) = H + 0 (x) + δH + 1 (x) + O(δ 2 ) , x ∈ Ω + (3.1) H - (δ) (x) = H - 0 (x; δ) + δH - 1 (x; δ) + O(δ 2 ) , x ∈ Ω - (3.2) with H - j (x; δ) = χ(y 3 ) H j (y β , y 3 δ ) . (3.3) In (3.1)-(3.
2), the symbol O(δ 2 ) means that the remainder is uniformly bounded by δ 2 , and in (3.3), the function y → χ(y 3 ) is a smooth cut-off with support in U -and equal to 1 in a smaller tubular neighborhood of Σ. The vector fields H j : (y α , Y 3 ) → H j (y α , Y 3 ) are profiles defined on Σ × R + : They are exponentially decreasing with respect to Y 3 and are smooth in all variables.

3.1. First terms of asymptotics in the conductor region. Hereafter, we present the construction of the first profiles H j = (H α j , h j ) and of the first terms H + j . The normal component h 0 of the first profile in the conductor is zero:

h 0 = 0 .
Then, the first term of the magnetic field in the dielectric region solves Maxwell equations with perfectly conducting conditions on Σ:

(3.4)      curl curl H + 0 -κ 2 H + 0 = curl j in Ω + H + 0 • n = 0 and curl H + 0 × n = 0 on Σ H + 0 × n = 0 and div H + 0 = 0 on Γ.
Thus the trace h 0 of H + 0 on the interface Σ is tangential. The first profile in the conductor region is exponential with the complex rate λ such that λ 2 = -iκ 2 :

(3.5) H 0 (y β , Y 3 ) = h 0 (y β ) e -λY 3 with λ = κ e -iπ/4 .
Note that, if h 0 is not identically 0, there exists C 0 > 1 independent of δ such that

(3.6) C -1 0 √ δ H - 0 ( • ; δ) 0,Ω -C 0 √ δ .
The next term which is determined in the asymptotics is the normal component h 1 of the profile H 1 :

(3.7) h 1 (y β , Y 3 ) = λ -1 D α h α 0 (y β ) e -λY 3 .
Here D α is the covariant derivative on Σ and we use the summation convention of repeated indices. The next term in the dielectric region solves:

(3.8)      curl curl H + 1 -κ 2 H + 1 = 0 in Ω + H + 1 • n = h 1 and curl H + 1 × n = iλh 0 on Σ H + 1 × n = 0 and div H + 1 = 0 on Γ.
Like above, h 1 is the trace of H + 1 on the interface Σ, and h α 1 denote its tangential components. The tangential components H α 1 of the profile H 1 are given by (3.9)

H α 1 (y β , Y 3 ) = h α 1 + Y 3 H h α 0 + b α σ h σ 0 (y β ) e -λY 3 , α = 1, 2 .
Here b γ α = a γβ b βα , where a γβ is the inverse of the metric tensor a αβ in Σ, and b αβ is the curvature tensor in Σ and

H = 1 2 b α α
is the mean curvature of the surface Σ. In particular, the sign of H depends on the orientation of the surface Σ. As a convention, the unit normal vector n on the surface Σ is inwardly oriented to Ω -, see Figure 2.

3.2. Asymptotic behavior of the skin depth. In a one-dimensional model, when the conductor Ω -is a half-space, the classical skin depth parameter is given by

(3.10) (σ) = 2 ωµ 0 σ .
This length corresponds to the distance from the surface of the conductor where the field has decreased of a rate e. In our situation, following [4, def. 4.1], we extend this definition to curved interfaces. For a data j, let us define

H (δ) (y α , y 3 ) := H - (δ) (x), y α ∈ Σ, 0 y 3 < h 0 , for h 0 small enough. Hereafter for any Z = (z 1 , z 2 , z 3 ) ∈ C 3 , |Z| denotes the vector- norm (|z 1 | 2 +|z 2 | 2 +|z 3 | 2 ) 1/2 in C 3 and • ,
• the corresponding hermitian scalar product. Definition 3.1. Let Σ be a smooth surface, and j a data of problem (2.1). such that for all y α in Σ, H (δ) (y α , 0) = 0. The skin depth is the length L(σ, y α ) defined on Σ and taking the smallest positive value such that

(3.11) |H (δ) y α , L(σ, y α ) | = |H (δ) (y α , 0)| e -1 .
Thus the length L(σ, y α ) is the distance from the interface where the field has decreased of a fixed rate. It depends on the conductivity σ and of each point y α in the interface Σ. A priori it also depends on the data j.

As a consequence of (3.5) and (3.9), there holds (3.12)

         |H (δ) (y α , y 3 )| 2 = |h 0 (y α )| 2 m(y α , y 3 ; δ) e -2y 3 Re(λ)/δ , with m(y α , y 3 ; δ) := 1 + 2y 3 H(y α ) + 2δ Re h 0 (y α ), h 1 (y α ) |h 0 (y α )| 2 + O (δ + y 3 ) 2 .
Relying on this formula, one can exhibit the asymptotic behavior of the skin depth L(σ, y α ) for high conductivity σ, cf. [4, Th. 4.2]:

Theorem 3.2. Let Σ be a regular surface with mean curvature H. Recall that (σ) is defined by (3.10). We assume that h 0 (y α ) = 0. The skin depth L(σ, y α ) has the following behavior for high conductivity:

(3.13) L(σ, y α ) = (σ) 1 + H(y α ) (σ) + O(σ -1 ) , σ → ∞ .
Remark 3.3. The higher order terms O(σ -1 ) in equation (3.13) do depend on the data j of problem (2.1).

AXISYMMETRIC DOMAINS

In order to perform scalar two dimensional computations which could represent correctly the features of a three-dimensional problem, we choose to consider an axisymmetric configuration in which Ω + and Ω -are axisymmetric domains with the same axis Ξ 0 : in cylindrical coordinates (r, θ, z) associated with this axis, there exists bi-dimensional "meridian" domains Ω m and Ω m ± such that

Ω = {x ∈ R 3 | (r, z) ∈ Ω m , θ ∈ T}, Ω ± = {x ∈ R 3 | (r, z) ∈ Ω m ± , θ ∈ T}.
Here T = R/(2πZ) is the one-dimensional torus. We denote by Γ m and Σ m the meridian curves corresponding to Γ and Σ, respectively, and by Γ 0 , Γ + 0 the following subsets of the rotation axis Ξ 0 , see Fig.

3

Γ 0 = Ξ 0 ∩ Ω m and Γ + 0 = Ξ 0 ∩ Ω m + . O r z Ω m - Ω m + Γ + 0 Γ m Γ 0 Σ m FIGURE 3. The meridian domain Ω m = Ω m -∪ Ω m + ∪ Σ m with boundary ∂Ω m = Γ m ∪ Γ 0
On such an axisymmetric configuration we consider a modification of problem (2.1)-(2.2): We take j ≡ 0 and impose instead of (2.2) non-homogeneous magnetic boundary conditions

(4.1) E • n = 0 and H × n = G × n on Γ , for a given data G ∈ H(curl, Ω) := {u ∈ L 2 (Ω) | curl u ∈ L 2 (Ω)}.
Then, the variational problem for the magnetic field H (δ) solution of the Maxwell equations (2.1)-(4.1) writes

Find H (δ) ∈ H 0 (curl, Ω) + G , such that for all K ∈ H 0 (curl, Ω), (4.2
)

Ω 1 ε(δ) curl H (δ) • curl K -κ 2 H (δ) • K dx = 0 . 4.1. Formulation in cylindrical components. For a vector field H = (H 1 , H 2 , H 3 ), we introduce its cylindrical components (H r , H θ , H z ) according to    H r (r, θ, z) = H 1 (x) cos θ + H 2 (x) sin θ , H θ (r, θ, z) = -H 1 (x) sin θ + H 2 (x) cos θ , H z (r, θ, z) = H 3 (x) ,
and we set

H(r, θ, z) = H r (r, θ, z), H θ (r, θ, z), H z (r, θ, z) .
We say that H is axisymmetric if H does not depend on the angular variable θ.

The Maxwell problem (4.2) is axisymmetric, which means that, expressed in cylindrical variables (r, θ, z) and components H, its coefficients do not depend on θ [START_REF] Hiptmair | Computation of resonant modes for axisymmetric Maxwell cavities using hp-version edge finite elements[END_REF][START_REF] Bernardi | Spectral methods for axisymmetric domains[END_REF][START_REF] Nkemzi | On the solution of Maxwell's equations in axisymmetric domains with edges[END_REF]. Recall that for a vector field H = (H 1 , H 2 , H 3 ) the cylindrical components of its curl write

(4.3)      (curl H) r = 1 r ∂ θ H z -∂ z H θ , (curl H) θ = ∂ z H r -∂ r H z , (curl H) z = 1 r ∂ r (rH θ ) -∂ θ H r .
As a consequence, if the right-hand side G is axisymmetric, and if (4.2) has a unique solution, then this solution is axisymmetric. According to (4.3), when H is axisymmetric its curl reduces to (4.4)

     (curl H) r = -∂ z H θ , (curl H) θ = ∂ z H r -∂ r H z , (curl H) z = 1 r ∂ r (rH θ ) . 4.2.
Axisymmetric orthoradial problem. We say that H is orthoradial if its components H r and H z are zero.

We assume that G is axisymmetric and orthoradial, i.e., G(r, θ, z) = (0, g(r, z), 0).

The components (H r , H z ) and H θ being uncoupled in (4.4), and the solution of problem (4.2) being unique, we obtain that this solution is orthoradial

H(δ) (r, θ, z) = (0, h (δ) (r, z), 0).

Variational formulation.

In this framework, the change from Cartesian to cylindrical coordinates requires the modification of the solution spaces H 0 (curl, Ω) used in problem (4.2). Precisely, the weighted space characterizing the orthoradial component

h (δ) (r, z) is V 1 1,Γ m (Ω m ) = {v ∈ H 1 1 (Ω m ) | v ∈ L 2 -1 (Ω m ) and v = 0 on Γ m }. Here, H 1 1 (Ω m ) = {v ∈ L 2 1 (Ω m ) | ∂ j r ∂ 1-j z v ∈ L 2 1 (Ω m ), j = 0, 1},
and for all α ∈ R, the space L 2 α (Ω m ) is the set of measurable functions v(r, z) such that

v 2 L 2 α (Ω m ) = Ω m |v| 2 r α drdz < +∞ . Remark 4.1. The space V 1 1,Γ m (Ω m ) incorporates essential boundary conditions, in partic- ular on Γ 0 , where v = 0, see [1, Remark II.1.1].
This leads us to solve the following two-dimensional scalar problem set in Ω m .

Find h (δ) ∈ V 1 1,Γ m (Ω m ) + g such that for all w ∈ V 1 1,Γ m (Ω m ), (4.5) a δ (h (δ) , w) = 0 ,
where in the conducting part Ω -:

a δ (h, w) := Ω m 1 ε(δ) ∂ z h ∂ z w + 1 r ∂ r (rh) 1 r ∂ r (rw) rdrdz -κ 2 Ω m h w rdrdz . 4.2.2. Asymptotic expansion. Let (r(ξ), z(ξ)) = τ (ξ), ξ ∈ (0, L), be an arc-length coordinate on the interface Σ m . Here ξ → τ (ξ) is a C ∞ function,
h + (δ) (r, z) = h + 0 (r, z) + δh + 1 (r, z) + O(δ 2 ) , h - (δ) (r, z) = h - 0 (r, z; δ) + δh - 1 (r, z; δ) + O(δ 2 ), h - j (r, z; δ) = χ(y 3 ) h θ j (ξ, y 3 δ ) .
Here the profiles h θ j are defined on Σ m × R + . Hereafter, we focus on the first terms h + 0 , h θ 0 , h + 1 , and h θ 1 . We introduce the interior and boundary operators

(4.7) D(r, z; ∂ r , ∂ z ) = ∂ 2 r + 1 r ∂ r + ∂ 2 z - 1 r 2 and B(ξ; ∂ r , ∂ z ) = -z (ξ)(∂ r + 1 r ) + r (ξ)∂ z .
In appendix B.1, we give the expansion of these operators in power series of δ inside the domain Ω -and on the interface Σ m . The terms h + 0 , h θ 0 , h + 1 , and h θ 1 satisfy the following problems coupled by their boundary conditions on the interface Σ m (corresponding to Y = 0) -compare with (3.4)-(3.5), (3.8)-(3.9):

(4.8)    Dh + 0 + κ 2 h + 0 = 0 in Ω m + , Bh + 0 = 0 on Σ m , h + 0 = g on Γ m ∪ Γ + 0 , (4.9) (∂ 2 Y -λ 2 )h θ 0 = 0 for 0 < Y < +∞ , h θ 0 = h + 0 for Y = 0 ,
with λ defined by (3.5), from which we deduce:

(4.10) h θ 0 (ξ, Y ) = e -λY h + 0 τ (ξ) . The next problem in the dielectric part is (4.11)    Dh + 1 + κ 2 h + 1 = 0 in Ω m + , Bh + 1 = -i∂ Y h θ 0 on Σ m , h + 1 = 0 on Γ m ∪ Γ + 0 , Thus (4.10) yields that h + 1 satisfies Bh + 1 = iλh + 0 Σ m on the interface. The next problem in the conductor part is (4.12) (∂ 2 Y -λ 2 )h θ 1 = -A 1 h θ 0 for 0 < Y < +∞ , h θ 1 = h + 1 for Y = 0 . Here, A 1 h θ 0 = -k + z r (ξ)∂ Y h θ 0 . From (4.10), we infer A 1 h θ 0 = e -λY λ k + z r (ξ)h + 0 Σ m .
Then, from equation (4.12), we obtain (4.13)

h θ 1 (ξ, Y ) = e -λY h + 1 (τ (ξ)) + Y 2 k + z r (ξ) h + 0 (τ (ξ)) .
Note that we can also deduce this profile h θ 1 from equations (3.5) and (3.9), but this is not obvious because the cylindrical coordinates are not a normal coordinate system, see Appendix B.2. Remark 4.2. Subsequently, we assume that the data g is a real valued function. Thus, the right hand side of the boundary value problem (4.8) is real. Hence h + 0 is a real valued function. Recall that Bh + 1 = κ e iπ/4 h + 0 Σ . From the boundary value problem (4.11), we infer: Re h + 1 = Im h + 1 . We will exploit this relationship in the numerical simulations of skin effect, see §6.

Configurations chosen for computations. We consider three classes of geometric configurations: one cylindrical configuration (A) and two spheroidal configurations (B and C).

Configuration A: cylindrical geometry. We assume that Ω is a circular cylinder of radius r 1 and length 1 , and Ω -is a coaxial cylinder of radius r 0 and length 0 . Hence, Ω m is a rectangle of width r 1 and length 1 , and Ω m -is a coaxial rectangle of width r 0 and length 0 , see Figure 4. We choose the parameters r 0 = 1, 0 = 2, r 1 = 2, 1 = 4 in computations. Right hand sides of problems. In configurations A and B, we take the data g = r. Hence, h (δ) satisfies the following inhomogeneous Dirichlet boundary condition

h (δ) (r, z) = r on Γ m .
For configuration C1, we take g = 0 and an interior data f := (curl j) θ with support inside Ω m + : f = 10 2 if r 2 /4 + z 2 0.8 and f = 0 otherwise.

FINITE ELEMENT DISCRETIZATIONS AND COMPUTATIONS

5.1. Finite element method. In this section, we consider two benchmarks for the computational domain: configuration A, see Figure 4, and configuration B, see Figure 5.

We use high order elements available in the finite element library Mélina, see [START_REF] Mélina | bibliothèque de calculs éléments finis[END_REF], and quadrangular meshes in the meridian domain. We discretize the variational problem (4.5). The script used to solve the problem (4.5) is adapted from [START_REF] Boissoles | Problèmes mathématiques et numériques issus de l'imagerie par résonance magnétique nucléaire[END_REF]. In the computations, we fix the angular frequency ω = 3.10 7 . We denote by h p,M (δ) the computed solution of the discretized problem (4.5) with an interpolation degree p and a mesh M. We define

A p,M σ := h p,M (δ) L 2 1 (Ω m -)
with σ = ωε 0 δ -2 , cf. (2.5) .

Interpolation degree.

We first check the convergence when the interpolation degree of the finite elements increases.

Configuration B1.

We consider the discretized problem with different degrees:

Q p 2
, for all p = 1, • • • , 20, and with three different meshes M 1 , M 3 and M 6 with 1, 3 or 6 layers of elongated elements in the skin region of the conductor Ω m -, see Figure 7.

2 Recall that Q p is the vector space of polynomials of two variables and partial degree p defined on the We represent in Figure 8 the absolute value of the difference between the weighted norms A p,M 1 σ and A 20,M 6 σ , versus p in semilogarithmic coordinates, and in each case: σ = 5 with circles, σ = 20 with squares, and σ = 80 with diamonds. We refer to [START_REF] Schwab | The p and hp versions of the finite element method for problems with boundary layers[END_REF] for theoretical results of convergence for the p-version in presence of an exponentially decreasing boundary layer. In Figure 9 we use mesh M 3 instead of M 1 . is the computed solution for several values of σ, with a fixed mesh M 3 (Figure 7). We represent | hσ | and | Im hσ | in configuration B1, with an interpolation degree p = 16, see Figure Figures 17 and18 show that the skin depth is larger for a fixed conductivity when the mean curvature of the conducting body surface is larger. Moreover, the sign of the

reference element K := [0, 1] × [0, 1]

POSTPROCESSING

In this section, we perform numerical treatments from computations in configuration B1, see subsection 6.2, and configuration A, see subsection 6.1, in order to investigate whether solutions are exponentially decreasing inside the conductor and with which rate. Let us recall that the standard skin depth (σ) is given by (3.10). Here the function s(y α , σ) := 1 ln 10

1 (σ) -H(y α )
depends on the skin depth (σ) and the mean curvature H(y α ) at the point y α of the interface Σ. We note that the mean curvature of the surface Σ is constant when z = 0. Hence, we introduce hereafter the theoretical slope s(σ) := s(y α , σ) when z = 0 .

In configuration B1, the principal curvatures at a point of the interface Σ when z = 0 are the constants κ 1 = 2 and κ 2 = 1 2 . Hence, H(y α ) = 5 4 . We infer

s(σ) = 1 ln 10 1 (σ) - 5 4 .
The accuracy of the asymptotic expansion is tested by representing the relative error between numerical and theoretical slopes:

err(σ) := s(σ) -s(σ) s(σ) ,
see the table 1. In order to make clear wether the influence of the curvature is visible in computations, we also display the theoretical curvature ratio curv ratio(σ) := The relative error decreases more than a half when the conductivity is multiplied by 4. Anyway, this relative error is much smaller than curv ratio(σ), which exhibits numerically the influence of the curvature on the skin depth. We perform similar computations with the mesh M 6 represented in Figure 7, see Table 2, and obtain still better results.

Remark 7.1. A similar postprocessing along the r-axis was performed for configuration A, see [START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]Ch. 8,§8.3.2]. Relative errors err(σ) are still consistent with the expansion (7.1). Non-radial postprocessings was performed along the segment OL, when L is the point with a colatitude ϕ 0 in the interface Σ m , see Figure 20: 1 for the theoretical values (σ), s(σ), and curv ratio(σ))

when ϕ 0 = π/2 - arctan (1/2) (i.e. r L = √ 2, z L = 1/ √ 2),
O L 2 4 r z Σ m ϕ 0 FIGURE 20.
In the meridian domain in configuration B1, L is the point on Σ m with a colatitude ϕ 0 = π/2 -arctan (1/2) 7.2. Configuration A. The mesh of the computational domain is the mesh M 4 represented in Figure 11. We extract values of log 10 | hσ | in Ω m -along the diagonal axis r = z, see Figure 21: here we denote by

ρ := (1 -r) 2 + (1 -z) 2 ,
the distance to the corner point a with coordinates (r = 1, z = 1), cf. Figure 4.

When compared with Figure 19, we see that the curves do not exactly behave like lines, which means that the exponential decay is not obvious. In order to measure a possible exponential decay, we define the slopes si (σ) of extracted values log 10 | hσ (ρ)| by:

si (σ) := log 10 | hσ (r i , z i )| -log 10 | hσ (z i+1 , r i+1 )| ρ i+1 -ρ i .
Here, ρ i is the distance to the corner point a defined by

ρ i := (1 -r i ) 2 + (1 -z i ) 2
with (r i , z i ) the extraction points. We present in the Figure 22 Whereas in configuration B1, the slopes clearly converge to a positive limit value as ρ = y 3 tends to 0, in configuration A the slopes tend to 0, which means that, stricto sensu, there is no exponential decay near the corner. Nevertheless we notice that in a region which is further away from the corner, a sort of exponential convergence is restored. This phenomenon is due to the fact that the principal asymptotic contribution inside the conductor is a profile globally defined on an infinite sector S (of opening π 2 in the present case) solving, instead the 1D problem (4.9) the model Dirichlet problem

(7.2) (∂ 2 X + ∂ 2 Y )h θ 0 -λ 2 h θ 0 = 0 in S , h θ 0 = h + 0 (a) on ∂S .

CONCLUSION

Even though addressing axisymmetric configurations for which the Maxwell system can be reduced to one scalar equation, our numerical experiments are in significative accordance with our theoretical results concerning the decay of solutions inside the conductor and their structure in the skin layer. Our asymptotics provide an a priori knowledge on solutions, which can be used for the design of meshes in view of a good quality finite element approximation: The mesh should fit the boundary of the conductor and can be coarse far from its boundary inside the conductor -depending on the skin depth. In this perspective, it is interesting to compare with [START_REF] Pardo | Two-dimensional highaccuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method[END_REF] where an adaptive a posteriori approach based on the hp method has been used for an industrial axisymmetric problem in electromagnetism.

APPENDIX A. ELEMENTS OF PROOF FOR THE MULTISCALE EXPANSION

Subsequently, we assume Assumption 2.1 on ω and Assumption 2.3 on the domains. In this framework, Theorem 2.2 gives the existence of δ 0 such that for all δ δ 0 , the problem (2.1)-(2.2) has a unique solution (E (δ) , H (δ) ) which is denoted by (E + (δ) , H + (δ) ) in the dielectric part Ω + , and (E - (δ) , H - (δ) ) in the conducting part Ω -. Furthermore, we suppose that the right hand side j ∈ H 0 (div, Ω) is smooth and its support does not meet the conductor domain Ω -.

Recall that (y α , y 3 ) is a local normal coordinate system to the surface Σ in U -, see Figure 2. The function y → χ(y 3 ) is a smooth cut-off with support in U -and equal to 1 in a smaller tubular neighborhood of Σ.

Theorem A.1. Under the above assumptions, the solution (E (δ) , H (δ) ) possesses the asymptotic expansion (see subsection A.7 below for precise estimates):

E + (δ) (x) ≈ j 0 δ j E + j (x) and H + (δ) (x) ≈ j 0 δ j H + j (x) , (A.1) E - (δ) (x) ≈ j 0 δ j E - j (x; δ) with E - j (x; δ) = χ(y 3 ) E j (y β , y 3 δ ) , (A.2) H - (δ) (x) ≈ j 0 δ j H - j (x; δ) with H - j (x; δ) = χ(y 3 ) H j (y β , y 3 δ ) , (A.3)
where E j (y β , y 3 δ ) → 0 and H j (y β , y 3 δ ) → 0 when y 3 δ → ∞. Moreover, for any j ∈ N, there holds

(A.4) E + j , H + j ∈ H(curl, Ω + ) and E j , H j ∈ H(curl, Σ × R + ).
Hereafter, we present elements of proof of this theorem and details about the terms in asymptotics (A.2)-(A.3). In §A.1, we expand the "magnetic" Maxwell operators in power series of δ inside the boundary layer U -. We deduce in §A.2 the equations satisfied by the magnetic profiles, and derive explicitly the first ones in §A.3. Then in §A.4 and §A.5, we do the same for the electric profiles. As an alternative, we show how to deduce directly the magnetic profiles from the electrical ones in §A.6. In §A.7, we conclude to the validation of the asymptotic expansion with a convergence result.

A.1. Expansion of the operators. Integrating by parts in the magnetic variational formulation (2.7), we find the following Maxwell transmission problem

(A.5)              curl curl H + (δ) -κ 2 H + (δ) = curl j in Ω + curl curl H - (δ) -κ 2 (1 + i δ 2 )H - (δ) = 0 in Ω - curl H + (δ) × n = (1 + i δ 2 ) -1 curl H - (δ) × n on Σ H + (δ) × n = H - (δ) × n on Σ H + (δ) × n = 0 on Γ.
It is important to notice that, since κ = 0, it is a consequence of the above equations that div H (δ) = 0 in Ω.

Therefore, we have in particular the extra transmission condition

(A.6) H + (δ) • n = H - (δ)
• n on Σ. We denote by L(y α , h; D α , ∂ h

3 ) the 2d order Maxwell operator curl curl -κ 2 (1 + i δ 2 )I set in U -in a normal coordinate system. Here D α is the covariant derivative on the interface Σ, and ∂ h 3 is the partial derivative with respect to the normal coordinate y 3 = h. Let a αβ (h) be the metric tensor of the manifold Σ h , which is the surface contained in Ω -at a distance h of Σ, see Figure 2. According to [START_REF] Dauge | Comportement asymptotique à haute conductivité de l'épaisseur de peau en électromagnétisme[END_REF][START_REF] Faou | Elasticity on a thin shell: formal series solution[END_REF], the metric tensor in such a coordinate system writes (A.7) a αβ (h) = a αβ -2b αβ h + b γ α b γβ h 2 , and its inverse expands in power series of h

a αβ (h) = a αβ + 2b αβ h + O(h 2 ) .
With this metric, a three-dimensional vector field H can be split into its normal component h and its tangential component that can be alternatively viewed as a vector field H α or a one-form field H α with the relation

(A.8) H α = a αβ (h)H β .
Subsequently, we use a property of the covariant derivative, that it acts on functions like the partial derivative:

D α h = ∂ α h. We denote by T(y α , h; D α , ∂ h 3 ) the tangent trace operator (1 + i δ 2 ) -1 curl • × n on Σ. If H = (H α , h), then (A.9) T(y α , h; D α , ∂ h 3 )H = (1 + i δ 2 ) -1 (∂ h 3 H α -D α h)dy α ,
From equations (iv)-(vi) in (A.32) for n = 0, and from (A.34), the first asymptotic of the electric field in the dielectric part solves the following problem

   curl curl E + 0 -κ 2 E + 0 = iωµ 0 j in Ω + E + 0 × n = 0 on Σ E + 0 • n = 0 and curl E + 0 × n = 0 on Γ.
According to the spectral Assumption 2.1, this boundary value problem is well-posed.

The next term determined in the asymptotic expansion is

E 1 . From equations (i)-(ii) in (A.32) for n = 1, E 1 satisfies for Y 3 ∈ I ∂ 2 3 E 1 (., Y 3 ) -λ 2 E 1 (., Y 3 ) = 0 ∂ 3 E 1 (., 0) = (curl E + 0 × n)(. , 0) . We denote by j k (y β ) = λ -1 (curl E + k × n)(y β , 0) for k = 0, 1. Hence, (A.35) E 1 (y β , Y 3 ) = -j 0 (y β ) e -λY 3 .
From equation (iii) in (A.32) for n = 1, and from (A.34), we obtain e 1 = 0. From equations (iv)-(vi) in (A.32) for n = 1, and from (A.35), the asymptotic of order 1 for the electric field in the dielectric part solves :

   curl curl E + 1 -κ 2 E + 1 = 0 in Ω + E + 1 × n = -j 0 × n on Σ E + 1 • n = 0 and curl E + 1 × n = 0 on Γ.
Then, from equations (i)-(ii) in (A.32) for n = 2, E 2,α solves the ODE for Y 3 ∈ I:

∂ 2 3 E 2,α (. , Y 3 ) -λ 2 E 2,α (. , Y 3 ) = -2b σ α ∂ 3 E 1,σ (. , Y 3 ) + b β β ∂ 3 E 1,α (. , Y 3 ) ∂ 3 E 2,α (. , 0) = (curl E + 1 × n) α (. , 0
) . We denote by j k,α the surface components of j k , for k = 0, 1. We obtain

E 2,α (y β , Y 3 ) = -j 1,α + λ -1 + Y 3 b σ α j 0,σ -H j 0,α (y β ) e -λY 3 . (A.36)
From equation (iii) in (A.32) for n = 2, we infer e 2 (y β , Y 3 ) = -λ -1 D α j α 0 (y β ) e -λY 3 . A. [START_REF] Faou | Elasticity on a thin shell: formal series solution[END_REF]. From E to H in the conducting part. An alternative way of calculating the magnetic profiles is to deduce them from the electric ones by means of a normal parameterization of the intrinsic curl operator, see [START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]Ch. 3]. Let E = (E α , e) be a 1-form fields in U -. There holds:

(A.37) (∇ × E) α = 3βα (∂ h 3 E β -∂ β e) and (∇ × E) 3 = 3αβ D h α E β on Σ h . Here, D h
α is the covariant derivative on Σ h , and is the Levi-Civita tensor, see [START_REF] Levi-Civita | The absolute differential calculus[END_REF][START_REF] Goldstein | Classical mechanics[END_REF]. The contravariant components ijk of depend on the normal coordinate h, and write in a normal coordinate system ijk = det a αβ (h) -1/2 0 (i, j, k) .

APPENDIX B. THE MULTISCALE EXPANSION OF THE ORTHORADIAL COMPONENT

We denote by ( e r , e θ , e z ) the basis associated with the cylindric coordinates (r, θ, z). B.1. Expansion of the operators. In the basis ( e r , e z ), recall that r(ξ), z(ξ) = τ (ξ), ξ ∈ (0, L) is an arc-length coordinate on the interface Σ m , and (ξ, h = y 3 ) is the associate normal coordinate system, see §4.2.2. The normal vector n(ξ) at the point τ (ξ) writes n(ξ) = -z (ξ), r (ξ) .

Hence, the tubular neighborhood U m -of Σ m inside Ω m -is represented thanks to the parameterization Ψ : (ξ, h) -→ (r, z) , (B.1) where r = r(ξ) -hz (ξ) , and z = z(ξ) + hr (ξ) .

Remark B.1. The curvature k(ξ) at the point τ (ξ) is defined by To perform the formal expansion of h - δ solution of (4.5) in Ω m -, we first use the change of variables Ψ in order to write the equations in the cylinder T L × [0, h 0 ). We then perform the rescaling

(B.2) k(ξ) = (r z -z r )(ξ) . For h 0 < 1/ k ∞ the change of coordinates Ψ is a C ∞ -diffeomorphism from the cylinder T L × [0, h 0 ) into U m -: U m -= Ψ T L × [0,
(B.3) Y = δ -1 h
in the equations set in Ω m -and Σ m in order to make appear the small parameter δ in the equations. Actually δ appears in the equations set in Ω -through the expression of the operator (1 + i δ 2 ) -1 D, where D is defined by (4.7). We obtain the formal expansion:

D = δ -2 ∂ 2 Y + δD 1 + δ 2
B = δ -1 ∂ Y - z r (ξ)
on the interface Σ m , and

(1 + i δ 2 ) -1 B = -iδ∂ Y + δ 2 P δ . Note that these expansions correspond to the expansions (A.9) in orthoradial symmetry. B.2. From the 3D asymptotic expansion to the 1D expansion. In this subsection, we show that we can obtain the profiles h θ 0 , and h θ 1 defined by (4.10) and (4.13) from the general profiles H α 0 and H α 1 defined above. We set (y α ) := (ξ, θ) a coordinate system on Σ. Thus, the normal coordinate system (y α , h) on U -is induced by the normal coordinate system (ξ, h) on U m -. The tubular neighborhood U -is parameterized by Φ : (y α , h) -→ (r cos θ, r sin θ, z) , where r and z are defined by (B.1). The associated tangent coordinate vector fields are The normal coordinate vector field is x 3 (h) = n(ξ) = -z (ξ) e r + r (ξ) e z .

A vector field H : (y α , h) → H(y α , h) in U -writes (here α = 1, 2) where k = k(ξ) is defined by (B.2). The inverse of the metric tensor a αβ in Σ is diagonal and there holds a 11 = 1 and a 22 = r -2 . From (B.10), we deduce the lemma.

H = H α x α (h) + h θ n ,

  and L is the length of the curve Σ m . Let (ξ, y 3 ) be the associate normal coordinate system in a tubular neighborhood of Σ m inside Ω m -. Then the normal vector n(ξ) at the point τ (ξ) can be written as (Frenet frame) (4.6) n(ξ) = -z (ξ), r (ξ) with z (ξ) = dz dξ and r (ξ) = dr dξ . Finally we denote by k(ξ) the curvature of Σ m in τ (ξ). In accordance with (3.1)-(3.2)-(3.3), we can exhibit series expansions in powers of δ for the magnetic field h (δ) which we denote by h + (δ) in the dielectric part Ω + , and by h - (δ)

FIGURE 4 .FIGURE 5 .FIGURE 6 .

 456 FIGURE 4. The meridian domain Ω m in configuration A

FIGURE 7 .

 7 FIGURE 7. The meshes M 1 , M 3 and M 6 for configuration B1

FIGURE 8 .

 8 FIGURE 8. Graph of A p,M 1 σ

FIGURE 9 .

 9 FIGURE 9. Graph of A p,M 3 σ

FIGURE 12 . 2 σ-A 16 ,M 3 σ

 122163 FIGURE 12. Graph of A p,M 2 σ

FIGURE 13 . 6 . 1 .

 1361 FIGURE 13. Graph of A 2,M k σ

  .

6. 2 .

 2 Skin effect in configuration B. Here hσ = h p,M 3 (δ)

FIGURE 14 .

 14 FIGURE 14. Configuration A. On the left, | hσ | when σ = 5. On the right, | hσ | when σ = 80

FIGURE 16 .

 16 FIGURE 16. Configuration B1. At the top, | hσ | when σ = 5 (on the left), and σ = 80 (on the right). At the bottom, | Im hσ | when σ = 5 and σ = 80

FIGURE 17 . 5 FIGURE 18 .

 17518 FIGURE 17. On the left, | Im hσ | in configuration B1 (H > 0). On the right, | Im hσ | in configuration C1 (H < 0). σ = 5

7. 1 .

 1 Configuration B. The mesh of the computational domain is the mesh M 3 , see Figure 7. We extract values of | hσ | in Ω m -along edges of the mesh M 3 for z = 0: in this configuration, the normal coordinate writes y 3 := 2 -r. Then, we perform a linear regression from values of log 10 | hσ (y 3 )| in the skin depth (σ), see Figure 19. We denote by n(σ) the number of extracted values on the axis z = 0 in the skin depth (σ). From the linear regression, we derive a numerical slope s(σ) such that log 10 | hσ (y 3 )| = -s(σ)y 3 + b , b ∈ R .

FIGURE 19 .

 19 FIGURE 19. In circles, extracted values of log 10 | hσ (y 3 )|. In solid line, linear regression of log 10 | hσ (y 3 )|. On the left σ = 20 (interpolation degree p = 12). On the right, σ = 80 (degree p = 16). n(20) = 6 , s(20) = 7.86906 and n(80) = 5 , s(80) = 16.29634. Configuration B1

  the graph of the slopes si (σ) for each curve on the Figure 21 (σ = 20, 80). For the sake of comparison, we also represent on the same figure the slopes of Figure 19 corresponding to configuration B1.

FIGURE 21 .

 21 FIGURE 21. In circles, extracted values of log 10 | hσ (ρ)|, ρ is the distance to the corner. On the left σ = 20 (interpolation degree p = 12). On the right, σ = 80 (degree p = 16). Configuration A

  h 0 ) . Thanks to (B.1), we obtain∂ r ∂ z = (1 -hk(ξ)) -1 r -(z + hr ) z (r -hz ) ∂ ξ ∂ h .

  R δ , where D 1 (ξ, Y ; ∂ ξ , ∂ Y ) = -k(ξ) + z r (ξ) ∂ Y ,and R δ is an operator, which has smooth coefficients in ξ and Y , bounded in δ. Hence, we obtaini(1 + i δ 2 ) -1 D + κ 2 I = A 0 + δA 1 + δ 2 Q δ , with A 0 = ∂ 2 Y -λ 2 I and A 1 = -k + z r (ξ) ∂ Y .Similarly, there holds

x 1

 1 (h) = (r (ξ) -hz (ξ)) e r + (z (ξ) + hr (ξ)) e z , (B.4) x 2 (h) = (r(ξ) -hz (ξ)) e θ . (B.5)

8 )+ b 2 σ h σ 0 - z r h 2 0eB. 3 .

 823 and corresponds to the magnetic field in cylindric components H(r, θ, z) = H r (r, θ, z) e r + H θ (r, θ, z) e θ + H z (r, θ, z) e z .Using (B.5) we easily obtain(r(ξ) -hz (ξ)) H 2 = H θ ,and introducing the stretched variable Y = δ -1 h, it follows(B.6) (r(ξ) -δY z (ξ)) H 2 = H θ .We insert the following expansionsH 2 = H 2 0 (y α , Y ) + δH 2 1 (y α , Y ) + O(δ 2 ), H θ = h θ 0 (ξ, Y ) + δh θ 1 (ξ, Y ) + O(δ 2 ) in equation (B.6). Then we perform the identification of terms with the same power in δ. We obtain the equations rAccording to (A.22), and (3.9) for α = 2, there holds(B.9) r(ξ)H 2 1 -Y z (ξ)H 2 0 = r(ξ) h 2 1 + Y H h 2 0 -λY .Lemma B.2. The main curvatures κ 1 , κ 2 of the interface Σ, and its mean curvatureH = 1 2 (κ 1 + κ 2 ) are κ 1 = k , κ 2 =We set r(ξ)h 2 j (y β ) = h + j (τ (ξ)), for j = 0, 1. Then, we infer (4.10)-(4.13) from the equations (A.22)-(3.9). Proof of the lemma B.2. The main curvatures are defined byκ 1 = b 1 1 , κ 2 = b 2 2, where b αβ is the curvature tensor. In this coordinate system on U -, the metric is defined byg ij (h) = x i (h), x j (h) R 3 . From (B.4)-(B.5), we obtain g ij (h) =   (r (ξ) -hz (ξ)) 2 + (z (ξ) + hr (ξ)) g αβ (h) = a αβ -2b αβ h + b γ α b γβ h 2 .Hence, the curvature tensor on the interface Σ is diagonal, and its diagonal components write (B.10) b 11 = k and b 22 = rz ,

  • 0,O the norm in L 2 (O). We quote from [3, Th.2.3]: Theorem 2.2. If the interface Σ is Lipschitz, under Assumption 2.1, there are constants σ 0 and C > 0, such that for all σ σ 0 , the Maxwell problem (2.1) with boundary condition (2.2) and data j ∈ H 0 (div, Ω) has a unique solution (E, H) in L 2 (Ω) 2 , which satisfies:

TABLE 1 .

 1 Postprocessing in configuration B1 with the mesh M 3

	5
	4 (σ) -5 1 4	.

TABLE 2 .

 2 see[START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF] Ch. 8, §8.5.1]. Relative errors are again very small. Postprocessing in configuration B1 with the mesh M 6 (see Table

	σ	5	20	80
	degree p	8	12	16
	n(σ)	13	9	7
	s(σ)	3.64239 7.88170 16.33051
	err(σ)	0.0084	0.001	0.0005

see [START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]Ch. 6, eq. ( 6.87)]. The operators L and T expand in power series of h with intrinsic coefficients with respect to Σ, see [START_REF] Faou | Elasticity on a thin shell: formal series solution[END_REF] for the 3D elasticity operator on a thin shell. We make the scaling Y 3 = δ -1 h to describe the boundary layer with respect to δ. Then, the three-dimensional harmonic Maxwell operators in U -are written L[δ] and T[δ]. These operators expand in power series of δ with coefficients intrinsic operators :

We denote by L n α and T n α the surface components of L n and T n . With the summation convention of repeated two dimensional indices (represented by greek letters), there holds

Here, ∂ 3 is the partial derivative with respect to Y 3 . We denote by L n 3 the transverse components of L n . There holds (A.12)

is the change of metric tensor. A.2. Equations for the coefficients of the magnetic field. The profiles H j and the coefficients H + j of the magnetic field satisfy the following system [START_REF] Naghdi | Foundations of elastic shell theory[END_REF] where I = (0, +∞). We perform in (A.13)-(A.14) the identification of terms with the same power in δ. The components of equation (A.13) are the collections of equations

for all n 2. Similarly, the surface components of the equation (A.14) write (A. [START_REF] Pardo | Two-dimensional highaccuracy simulation of resistivity logging-while-drilling (LWD) measurements using a self-adaptive goal-oriented hp finite element method[END_REF]) curl H + 0 × n = 0 , and

for all n 1. Using the expression of the operator L 0 , and expanding H + (δ) in Ω + , we thus see that, according to the system (A.5), the profiles H n = (H n , h n ) and the terms H + n have to satisfy, for all n ≥ 0, (A.17)

where λ = κ e -iπ/4 , cf (3.5) (so that -λ 2 = iκ 2 ) and h n denotes the trace of H + n on Σ. In (A.17), we use the convention that the sums are 0 when n = 0. The transmission condition (A.6) implies the extra continuity condition

The set of equations (A.17)-(A.18) allows to determine H n and H + n by induction.

A.3. First terms of the magnetic field asymptotics. According to equation (i) in (A.17), the normal component h 0 of the first profile in the conductor vanishes:

(A. [START_REF] Schwab | The p and hp versions of the finite element method for problems with boundary layers[END_REF])

Hence, according to (A.17) (ii)-(iv) and (A.18), the first term of the magnetic field in the dielectric region solves Maxwell equations with perfectly conducting conditions on Σ:

Thus the trace h 0 of H + 0 on the interface Σ is tangential. According to equations (v)-(vi) in (A.17), H 0 satisfies the following ODE

The unique solution of (A.21) such that H 0 → 0 when Y 3 → ∞, is, with the choice (3.5) for λ, the tangential field H 0 (y β , Y 3 ) = h 0 (y β ) e -λY 3 . Combining with (A. [START_REF] Schwab | The p and hp versions of the finite element method for problems with boundary layers[END_REF], we find that the first profile in the conductor region is exponential with the complex rate λ:

The next term which is determined in the asymptotics is the normal component h 1 of the profile H 1 given by equation (i) of (A.17) for n = 1. We obtain

According to (A.17) (ii)-(iv) and (A.18), the next term in the dielectric region solves:

Recall that h 1 is the trace of H + 1 on the interface Σ. We denote by h 1,α its tangential components. According to equations (v)-(vi) in (A.17) for n = 1, H 1 satisfies the following ODE (for Y 3 ∈ I)

From (A.22), the unique solution of (A.25) such that

Using the relation (A.8) and performing the scaling h = δY 3 in the previous equation, we obtain for the contravariant components

), and (A.26), the tangential components H α 1 are given by (3.9)

Remark A.2. Note that the boundary value problems (A.20) and (A.24) are well-posed.

It is a consequence of the spectral Assumption 2.1 on ω.

A.4. Equations for the coefficients of the electric field. The second order Maxwell operator for the electric field writes

We denote by B(y α , h; D α , ∂ h 3 ) the tangent trace operator curl • × n on Σ in a normal coordinate system. If E = (E α , e), then

see [START_REF] Éron | Modélisation mathématique de phénomènes électromagnétiques dans des matériaux à fort contraste[END_REF]Ch. 3,Prop. 3.36]. We define B[δ] the operator obtained from B in U -after the scaling Y 3 = δ -1 h. This operator expands in power of δ :

Recall that ∂ 3 is the partial derivative with respect to Y 3 . Thus, denoting by B n α the surface components of B n , we obtain

According to the second and third equations in system (A.27), the profiles E j and the terms E + j of the electric field satisfy the following system

We identify in (A.29)-(A.30) the terms with the same power in δ. The components of the equation (A.29) are collections of equations, similar to the equations (A.15) set for the magnetic field. The surface components of the equation (A.30) write

According to the system (A.27) and (A.28), the profiles E n = (E n , e n ) and the terms E + n have to satisfy, for all n 0, (we recall λ = κ e -iπ/4 , cf (3.5), so that -

Hereafter, we determine the terms E n and E + n by induction. A.5. First terms of the asymptotics for the electric field. According to equations (i)-(ii) in system (A.32) for n = 0, E 0 satisfies the following ODE (A.33)

The unique solution of (A.33) such that E 0 → 0 when Y 3 → ∞, is E 0 = 0. From equation (iii) in system (A.32) for n = 0, there holds e 0 = 0. We infer

Here, a αβ (h) is the metric tensor of the manifold Σ h , see (A.7). The indices i, j, k ∈ {1, 2, 3}, and 0 (i, j, k) equals 1 when (i, j, k) is an even circular permutation, and equals -1 when (i, j, k) is an odd circular permutation, and 0 (i, j, k) = 0 otherwise.

Remark A.3. Let a = det(a αβ ), and recall that b ν ν = 2H. Then using (A.7) we obtain ijk = a -1/2 1 + 2Hh + O(h 2 ) 0 (i, j, k) .

We make the scaling h = δY 3 and expand equations (A.37) in power series of δ:

) . Here, j αβ = a -1/2 0 (α, β, 3) . From the expansions (A.2) and (A.1) combined with Faraday's law (2.1) written in normal coordinates, we obtain the profiles H j in the expansion (A.3) of the magnetic field from the profiles E j of the electrical field. In particular, from (A.34), (A.35), and (A.36), we obtain explicitly the first terms H 0 and H 1 , see (A.22), (3.9), and (A.23). A.7. Convergence result. The validation of the asymptotic expansion (A.1), (A.2),(A.3), and (A.4), consists in proving estimates for remainders defined as

By construction of the terms E j = (E + j , E - j ) in the dielectric and conductor parts, the remainders R E m; δ satisfy the assumption of [3, Th. 5.1], which is an estimate on the right hand side when the Maxwell operator is applied to R E m; δ . Thus [START_REF] Caloz | Uniform estimates for transmission problems with high contrast in heat conduction and electromagnetism[END_REF]Th. 5.1] yields that for all m ∈ N, and δ ∈ (0, δ 0 ), there holds the optimal estimate (A.41) R E,+ m; δ 0,Ω + + curl R E,+ m; δ 0,Ω + + δ -1 2 R E,- m; δ 0,Ω -+ δ In Ω -, the presence of profiles in the expansion prevents to link R H,- m; δ and R E,- m; δ via the Maxwell equations in a similar way as (A.42). Nevertheless, there holds also uniform estimates for R H,- m; δ : δ -1 2 R H,- m; δ 0,Ω -+ δ


Paolo A Calligari
email: calligari@cnrs-orleans.fr

Gerald R Kneller

Andrea Giansanti

Paolo Ascenzi

Alessandro Porrello

Alessio Bocedi
email: bocedi@uniroma3.it

Inhibition of viral group-1 and group-2 neuraminidases by oseltamivir: a comparative structural analysis by the ScrewFit algorithm

Keywords: flu, influenza, H, hemagglutinin, NA, neuraminidase neuraminidase, oseltamivir binding, enzyme structure, ScrewFit algorithm, comparative structural analysis

The viral surface glycoprotein neuraminidase (NA) allows the influenza virus penetration and the egress of virions. NAs are classified as A, B, and C. Type-A NAs from influenza virus are subdivided into two phylogenetically distinct families, group-1 and group-2. NA inhibition by oseltamivir represents a therapeutic approach against the avian influenza virus H5N1. Here, structural bases for oseltamivir recognition by group-1 NA1, NA8 and group-2 NA9 are highlighted by the ScrewFit algorithm for quantitative structure comparison.

Oseltamivir binding to NA1 and NA8 affects the geometry of Glu119 and of regions Arg130-Ser160, Val240-Gly260, and Asp330-Glu382, leading to multiple NA conformations.

Additionally, although NA1 and NA9 share almost the same oseltamivir-bound final conformation, they show some relevant differences as suggested by the ScrewFit algorithm.

These results indicate that the design of new NA inhibitors should take into account these family-specific effects induced on the whole structure of NAs.

A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T

http://www.who.int/csr/disease/influenza/en/). Flu viruses, belonging to the Orthomyxoviridae family, are classified as A, B, and C, on the basis of the antigenic differences of their nuclear and matrix proteins,; all avian flu viruses belong to type-A. The antigenicity of the two surface glycoproteins hemagglutinin (H) and neuraminidase (NA) is used to cluster type-A flu viruses into subtypes: 16 for H (H1-H16) and 9 for NA (NA1-NA9). Virus uptake to the cell-surface sialic acid receptor is mediated by H. After virus replication, NA removes sialic acid from viral and cellular glycoproteins to facilitate virus release. Different combinations of H and NA subtypes (i.e., HxNy) have been found in avian species [START_REF] Horimoto | Pandemic threat posed by avian influenza A viruses[END_REF][START_REF] Webster | Evolution and ecology of influenza A viruses[END_REF][START_REF] Fouchier | Characterization of a novel influenza A virus hemagglutinin subtype (H16) obtained from black-headed gulls[END_REF].

The possible transmission of mutant viruses of avian flu to humans indicate that the danger of a pandemic is real, thus urging the scientific community to develop new anti-flu drugs, able to respond to the repertoire of antigenic mutations in the H5N1 virus strains [START_REF] Horimoto | Pandemic threat posed by avian influenza A viruses[END_REF][START_REF] Webster | Evolution and ecology of influenza A viruses[END_REF][START_REF] Fouchier | Characterization of a novel influenza A virus hemagglutinin subtype (H16) obtained from black-headed gulls[END_REF].

Two classes of antiviral drugs are available against flu viruses: the inhibitors of the ion channel activity of the M2 membrane protein (i.e., amantadine and rimantadine), and the NA inhibitors (i.e., oseltamivir (Tamiflu ® ) and zanamivir (Relenza ® )). The resistance against both oseltamivir and zanamivir reflects single point mutations of highly conserved NA residues [START_REF] Clercq | Antiviral agents active against influenza A viruses[END_REF][START_REF] Hurt | Susceptibility of highly pathogenic A (H5N1) avian influenza viruses to the neuraminidase inhibitors and adamantanes[END_REF]. Specifically, the sensitivity of flu viruses to oseltamivir, expressed by inhibition of the NA activity (IC 50 ), is different between oseltamivir-sensitive viruses (IC 50 ≅ 0.1-10 nM), and drug-resistant viruses (IC 50 ≅ 90 nM) [START_REF] Mai Le | Isolation of drug-resistant H5N1 virus[END_REF].

NAs are homo-tetrameric glycoproteins of ~ 240 kDa belonging to the three different types A, B, and C. Type-A NAs are subdivided in two phylogenetically distinct groups. Group-1 contains NA1, NA4, NA5, and NA8, whereas group-2 comprises NA2, NA3, NA6, NA7, and NA9 [START_REF] Kim | Influenza neuraminidase inhibitors possessing a novel hydrophobic interaction in the enzyme active site. Design, synthesis, and structural analysis of carbocyclic sialic acid analogues with potent anti-influenza activity[END_REF]. Notably, NA2 and NA9 have been used to design selective NA inhibitors oseltamivir and zanamivir [START_REF] Kim | Influenza neuraminidase inhibitors possessing a novel hydrophobic interaction in the enzyme active site. Design, synthesis, and structural analysis of carbocyclic sialic acid analogues with potent anti-influenza activity[END_REF].

Here, a detailed comparative analysis of the three-dimensional structures of

A C C E P T E D M A N U S C R I P T

oseltamivir-free and oseltamivir-bound viral NA1, NA8, and NA9 by the ScrewFit algorithm is reported. The latter represents a new quantitative approach for structure comparison based on the analysis of the screw motions of contiguous peptide planes along a protein backbone [START_REF] Kneller | Superposition of molecular structures using quaternions[END_REF][START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF][START_REF] Calligari | ScrewFit: a novel approach for continuum protein secondary structure assessments[END_REF].

Finally, we show that the ScrewFit analysis of NAs structures confirm most of the recent experimental results in respect to the differences between group-1 and group-2 NAs [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF]. In particular, we show that the ligand-free NA1 conformation is very similar to the oseltamivir-complexed NA1 open conformation whereas the closed ligand-bound conformation shows some relevant differences from them. We also display that these differences are group-dependent as they are also shown by the same comparison for NA8 conformations (group-1) but not in NA9 ones (group-2). This result is coherent to the fact that oseltamivir has been designed using group-2 NAs as templates [START_REF] Kim | Influenza neuraminidase inhibitors possessing a novel hydrophobic interaction in the enzyme active site. Design, synthesis, and structural analysis of carbocyclic sialic acid analogues with potent anti-influenza activity[END_REF].

Interestingly, the three-dimensional structure of the NA1 oseltamivir-bound closed conformation and NA9 are very similar, excepting some differences not reported in previous studies [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF] and detected here by means of ScrewFit. These results indicate that NA inhibitors should be designed expressly taking the target enzyme as the template; in this way crossinhibition of group-1 and group-2 NAs could not be excluded, resulting less efficacy by a effective therapeutically point of view.

Materials and Methods

Sequence analysis of NAs

The amino acid sequences of NAs from flu type-A viruses NA1 (PDB id: 2HTY), NA2 (PDB id: 2BAT), NA3 (NCBI id: AAO62020.1), NA4 (PDB id: 2HTV), NA5 (NCBI id: AAA43672.1), NA6 (PDB id: 1V0Z), NA7 (NCBI id: AAA43425.1), NA8 (PDB id: 2HT5), and NA9 (PDB id: 7NN9), and from the flu type-B virus (PDB id: 1INF) were recovered from the NCBI database (available at: http://www.ncbi.nlm.nih.gov), and from the Protein Data Bank [START_REF] Berman | The Protein Data Bank[END_REF]. The amino acid sequences have been multi-aligned with ClustalW (available at: http://www.ebi.ac.uk) (Fig. S1). The aminoacid residue numbering adopted here corresponds to that reported by Russell and coworkers [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF].

ScrewFit analysis of NA three-dimensional structures

The three-dimensional structures used in the ScrewFit analysis were obtained from the Protein Data Bank [START_REF] Berman | The Protein Data Bank[END_REF]: NA1 (ligand-free, PDB id: 2HTY; oseltamivir-bound open conformation, PDB id: 2HU0; oseltamivir-bound closed conformation, PDB id: 2HU4) [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF], NA8 (ligand-free, PDB id: 2HT5, oseltamivir-bound 'open' conformation, PDB id: 2HT7; oseltamivir-bound 'closed' conformation, PDB id: 2HT8) [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF], and NA9 (ligand-free, PDB id: 7NN9) [START_REF] Varghese | Three-dimensional structure of the complex of 4guanidino-Neu5Ac2en and influenza virus neuraminidase[END_REF]; oseltamivir-bound conformation (PDB id: 2QWK) [START_REF] Varghese | Drug design against a shifting target: a structural basis for resistance to inhibitors in a variant of influenza virus neuraminidase[END_REF]) (Fig. S2). The analysis is just limited to these structures for a better comprehension of the relationships among proteins belonging to type-A and -B NAs crystallized with the same compound (e.g., oseltamivir).

The ScrewFit algorithm [START_REF] Kneller | Superposition of molecular structures using quaternions[END_REF][START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF] superposes molecular structures, defined by atomic coordinates, combining quaternionic representation of rotation matrices and Chasles' theorem on rigid-body displacements. When applied to subsequent peptide planes in protein structures, ScrewFit gives local helical parameters of the protein backbone winding (see: http://dirac.cnrs-orleans.fr/Screwit/) [START_REF] Kneller | Superposition of molecular structures using quaternions[END_REF][START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF][START_REF] Calligari | ScrewFit: a novel approach for continuum protein secondary structure assessments[END_REF].

The ScrewFit algotithm takes into consideration three structural parameters: the orientational distance, the radius of screw motion, and the straightness. For each pair of consecutive peptide planes represented by their atoms {C, O, N}, ScrewFit defines three parameters indicating their relative

A C C E P T E D M A N U S C R I P T

orientation and distance from a common axis of rotation (the axis of screw motion) (Fig. 1).

The oxygen, carbon, and nitrogen atoms can be used to define the plane of a peptide group. Two subsequent peptide planes can be superposed as rigid bodies, through the minimization of the target function

m q ( ) = ω α D q ( ) ⋅ x α -x α ' ( ) 2 α=1 3 ∑ , (1)
where x α { } and ′

x α { } are the coordinates of the atoms in the reference and target plane, respectively. Both sets of coordinates are to be taken relative to the respective centers of rotation, which is here chosen to be the C-atom. The index α runs over the three atomic species {α = O, C, N} and to each atom is assigned a positive weight ω α value proportional to its mass (with ω α α=1 3 ∑ = 1 ). The vector q ≡ q 0 ,q 1 ,q 2 ,q 3 { } comprises four real normalized quaternion components, with q 0 2 + q 1 2 + q 2 2 + q 3 2 = 1 and D is a proper rotation matrix, with det(D)=+1, which has the form

D = q 0 2 + q 1 2 -q 2 2 -q 3 2
2(-q 0 q 3 + q 1 q 2 ) 2(q 0 q 2 + q 1 q 3 ) 2(q 0 q 3 + q 1 q 2 ) q 0 2 + q 2 2q 1 2q 3 2 2(-q 0 q 1 + q 2 q 3 ) 2(q 0 q 3 + q 1 q 3 ) 2(q 0 q 1 + q 2 q 3 ) q 0

2 + q 3 2 -q 1 2 -q 2 2 ⎛ ⎝ ⎜ ⎜ ⎜ ⎞ ⎠ ⎟ ⎟ ⎟ .
The target function ( 1) is to be minimized with respect to the four quaternion components. As demonstrated in previous work [START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF], this constrained minimization can be mapped onto an eigenvector problem of the form Mq = q, where M

≡ M x α ,x α ' { } (
) is a positive semi-definite matrix related to the target function. The solution of this kind of problem consists of four eigenvalues, λ j j = 1,K,4 ( ) , which coincide with the errors in the superposition fit described by the target function. The quaternion corresponding to the smallest eigenvalue, λ 1 = λ min , is the solution for the optimal fit which leads to the superposition of the two planes and, hence, it gives a description of the relative orientation of

′ x α { } with respect to x α { }. The largest A C C E P T E D M A N U S C R I P T
eigenvalue, λ 4 = λ max , describes, in contrast, the "worst" superposition and gives the maximal Euclidean distance between the two peptide planes. Through λ max , it is possible to define the first ScrewFit parameter, which is the unique orientational distance between two consecutive peptide planes,

Δ Ω = ω α x α -′ x α ( ) 2 α=1 3 ∑ λ max . ( 2
)
One notes that 0 ≤ Δ Ω ≤ 1.

The second parameter is radius of the screw motion,

ρ = | t ⊥ | 2 1+ cot 2 φ /2 ( ) . (3)
It is derived from Chasles' theorem, which states that any rigid-body displacement (like the superposition of two consecutive peptide planes) can be described by a rotation around an axis passing through a reference point, which may be located outside the rigid body, and a translation along the rotation axis. Here t ⊥ is the component of the distance vector t = CC', which is orthogonal to the rotation and translation axis n of the screw motion, and φ is the rotation angle.

Both n and φ are obtained through the relation

q = cos(φ /2),sin(φ /2)n x ,sin(φ /2)n y ,sin(φ /2)n z { } ,
where q is the quaternion describing the optimal rotation .

The straightness σ at residue i is defined by

σ i = μ i T ⋅ μ i+1 , (4)
where μ i is given by

μ i = R i+1 ⊥ -R i ⊥ | R i+1 ⊥ -R i ⊥ | . ( 5
)
The vector R i ⊥ in Eqn. ( 5) locates the point, on the helix axis, which is closest to the u ⊥ = 1/2(u ⊥ + cot(φ /2)n × t) . One recognizes that ρ = u ⊥ .

The radius ρ (see Eqn. Straightness σ is the scalar product between the unit vectors of the axis of screw motion, relative to four consecutive peptide planes, it gives information about local curvatures or kinks of a protein backbone tract.

It is worth to note here that both ρ and σ, defined as in ( 3) and ( 4

A C C E P T E D M A N U S C R I P T

Results

Oseltamivir binding to group-1 NAs

The comparative analysis of the ligand-free against the oseltamivir-bound open and closed conformations of NA1 indicates that there is not a remarkable difference between the ligand-bound open conformation and the ligand-free structure of NA1, while a significant difference is observed by comparison of the drug-free enzyme and the oseltamivir-bound closed conformation (Fig. 2). The changes find between the two NA1 oseltamivir-bound structures depend on different concentrations (20 μM and 0.5 mM, respectively) of oseltamivir and on different times of crystal soaking (30 min and 3 days, respectively) [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF]. As a whole, these observations indicate that oseltamivir binding induces a structural shrinking of NA1, mainly localized just in few regions.

Figure S3 report the ScrewFit profiles differences among NA8 ligand-free, open and closed conformations , showing that the ScrewFit parameters profiles are quite similar to the NA1 (Fig. 2).

The major structural changes, in both NA1 and NA8, are located in three regions:

Arg130-Ser160 (region A; also named 150-loop, see [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF]), Val240-Gly260 (region B), and Asp330-Glu385 (region C). The inspection of the aminoacid sequence multi-alignement indicates that these regions comprise highly conserved residues involved in ligand binding Glu, and Arg in region A, Ser, and Ala in region B, and Arg in region C (Fig. S1). Regions A and B are characterized by patterns of highly conserved residues, whereas in region C there is more variability. Patterns of mutations in these three regions differentiate group-1 and group-2 NAs, and are responsible for the different oseltamivir binding behavior (see below).

In addition to the regions mentioned above, small conformational change(s) are found also in a narrow region close to Glu119, which is located at the NA1 active site (Fig. 3A). which in group-1 NAs makes an hydrogen bond with His274 (Fig. 3C). His274 is known to be the site of the critical mutation His274Tyr, which leads to high resistance of group-1 NAs against oseltamivir [START_REF] Yen | Neuraminidase inhibitor-resistant influenza viruses may differ substantially in fitness and transmissibility[END_REF]. Interestingly, the open-to-closed transition induces a slight change in the length of the hydrogen bond between residues His274 and Pro245 (Fig. 3C). Moreover, the hydrogen bond between residues Gln249 and Ala270 is also affected by the open-toclosed conformational transition, shortening from 3.25 Å to 2.97 Å .

A C C E P T E

In contrast, oseltamivir binding does not alter significantly the direction of the screw motion axis in region C through the open-to-closed transition; indeed, the difference-plot of straightness parameter is quite flat in this region (Fig. 3D). Nevertheless, this region presents several variations in the orientation of the peptide planes as suggested by changes of the orientational distance and helix radius parameters (see Fig. 2). The latter undergo, between residue Asp330 and Pro340, similar variations as in region B. The rest of local changes in region C are dependent on structural modifications of the intra-chain hydrogen bonds, indeed, they are correlated to changes in the orientation of residues located in the outer protein surface [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF].

Oseltamivir binding to group-2 NAs

Here we will show the effects of oseltamivir binding to a representative member of group-2 NAs. Figure 4 shows the difference ScrewFit profile of native NA9 minus the

A C C E P T E D M A N U S C R I P T

enzyme:oseltamivir complex. Oseltamivir binding does not induce structural change(s) in this group-2 NA. This result agrees with the fact that oseltamivir was designed to optimally adapt to group-2 NAs [START_REF] Luo | Antiviral drugs fit for a purpose[END_REF].

Comparison between NA1 and NA9 binding with oseltamivir

In Fig. 5 are reported the difference ScrewFit plots of the open and closed configurations of NA1 minus the ligand-free structure of group-2 NA9. These difference plots are focalized on the regions where major structural changes have been found by the intragroup ScrewFit differences plot (Figures 2 and3). This analysis is suggested by the observation that the structure of group-1 NAs tends to resemble that of group-2 upon oseltamivir binding [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF]. Following this statement, it would be expected a zero difference The region close to the Glu119 seems to confirm a better superposition of oseltamivirbound NA1 closed conformation with NA9 free in respect to NA1 open conformation as observed in [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF] (Fig. 5A).

In Region A, differences are observed at sites Leu139-Leu140 and Ser145-Arg152, indicating that this region is influenced by oseltamivir binding to group-1 NAs. In particular, residues Gly147 and Arg152 change their configuration through the open-to-closed transition (Fig. 5B). Indeed, in the closed state, they assume a conformation similar to the one in NA9.

Nevertheless, in the same region some discrepancies between NA1 closed conformation and NA9 are still present. In particular residues Ser145 and Asp151 shows drastically different backbone conformations in the NAs. This local difference is also confirmed by the direct

A C C E P T E D M A N U S C R I P T

visual comparison of the structures shown in Fig. 5B. In Region B, there are pronounced deviations of all three ScrewFit parameters, at the site of the catalytically relevant Ser246.

This deviation almost disappears in the NA1 closed conformation (Fig. 5C). In Region C, although the pattern of structural differences is more confused, any significant change seems

to differentiate the open and the closed configurations of NA1 in respect to the ligand-free NA9 (Fig. 5D).

Discussion and Conclusions

The drug-dependent 'induced fit' transition of group-1 to group-2 NAs observed by Russell and coworkers [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF], is confirmed and refined here by using the ScrewFit algorithm developed for quantitative molecular structure comparison [START_REF] Kneller | Superposition of molecular structures using quaternions[END_REF][START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF][START_REF] Calligari | ScrewFit: a novel approach for continuum protein secondary structure assessments[END_REF]. Indeed, when crystals of NA1 are incubated with oseltamivir a transient open conformation different from that of the drug-free form and from that of the final stable oseltamivir-bound derivative is observed. In turn, the ligand free NA1 and NA9 structure are different, but the final stable oseltamivirbound NA1 and NA9 conformations are similar. When the oseltamivir-mediated open-toclosed transition in NA1 and NA8 takes place, both Glu119 and Asp151 move toward the bound oseltamivir (Fig. 3). Nevertheless, some significant differences between the oseltamivir-bound group-1 and group-2 are pinpointed by ScrewFit (Fig. 5) and refine the results reported in [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF].

Indeed, although the reported 'induced fit' is evident from Figures 2 and3, all regions present some sites with an opposite or neutral behavior in respect to this transition. In particular, the residue Asp151, which is known to interact with substrate, clearly assume different backbone configuration in the NA1 closed conformation and in NA9. These results

A C C E P T E D M A N U S C R I P T

show that the open-to-closed transition does not correspond to a homogeneous 'induced-fit' of group-1 to group-2 NAs [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF][START_REF] Luo | Antiviral drugs fit for a purpose[END_REF] and reveal a more complex adaptation of the whole threedimensional structure. Another important clue comes from Figure 4, where it is shown that NA9 structure is essentially unaffected by oseltamivir binding, suggesting that drug-binding has some thermodynamic optimality. It is worth noting that different studies were able to define the active site "geography" of NAs and the location of the residues that recognize the substrate [START_REF] Stoll | Influenza neuraminidase inhibitors: structure-based design of a novel inhibitor series[END_REF][START_REF] Wang | Comparative binding energy (combine) analysis of influenza neuraminidase-inhibitor complexes[END_REF][START_REF] Wang | Design, synthesis, and structural analysis of influenza neuraminidase inhibitors containing pyrrolidine cores[END_REF][START_REF] Re | Remarkable loop flexibility in avian influenza n1 and its implications for antiviral drug design[END_REF] (i.e., Arg118, Glu119, Asp151, Arg152, Trp178, Ile222, Arg224, Glu227, Ala/Ser246, Glu276, Glu277, Arg292 and Arg371) (Fig. S1, bold characters). Interestingly, most of the small kinks in the difference plots refer to these residues. These small structural changes could be further exploited in more refined analyses by the ScrewFit algorithm, especially in respect to the minor discrepancies found between NA1 and NA8.

Moreover, it is worth to note that all results from the direct NA1 and NA9 comparison (Fig. 5), have been reported here only qualitatively but they can also be used in a more quantitative manner as a measure of the local discrepancies between different conformations.

The difference plots of the open (induced by 20 μM oseltamivir) and closed (induced by 0.5 mM oseltamivir) configurations of NA1 minus ligand-free structure indicate that at mild inhibitor concentrations there are local distortions not only in the region of the 150-loop [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF], but also at residues Val240-Ala250, in the Asp330-Pro340 region and close to residues Pro380 and Phe410 (Fig. 2). Aminoacid mutations taking place in these regions, involved in ligand recognition, may modulate drug-binding.

These results confirm and complete the perspectives suggested in [START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF] and [START_REF] Luo | Antiviral drugs fit for a purpose[END_REF] about the design of new antiviral drugs in respect to the effects of inhibitors on the whole structure of NAs. The ScrewFit analysis, although related only on backbone conformational changes, proved to be useful in the evaluation and refinement of crystallographic results. Moreover it seems to be potentially relevant also in the assessment of protein structural homology, both locally

A C C E P T E D M A N U S C R I P T

and globally.

From a structural point of view, present results raise the basic question whether group-1

NAs are more rigid than group-2 enzymes in oseltamivir binding and which role the variability of the regions flanking the active site plays in respect to the flexibility of NAs.

Although answering to this kind of questions is not possible through a static analysis, and would require kinetic experiments and dynamic simulations, it is reasonable to assume that, through patterns of function conserving mutations, viral proteins explore not only structural possibilities, but also dynamical repertoires. This hypothesis seems to be confirmed by a preliminary molecular simulation made on ligand-free and oseltamivir-bound NA1 [START_REF] Re | Remarkable loop flexibility in avian influenza n1 and its implications for antiviral drug design[END_REF],

showing that NA1 could explore an ensemble of conformations, same of them even more squashed than the closed or wider than the open one. Therefore, important clues could emerge from the ScrewFit analysis of the NA4 and NA5 structures ligand-free and oseltasmivir-bound, when available. Additionally, it would be interesting to have further comparative molecular dynamics simulations of NA1 (or NA8) and NA9, with the aim of checking if phylogenetical distances correspond to different dynamical behaviors.

A C C E P T E D M A N U S C R I P T

  [START_REF] Fouchier | Characterization of a novel influenza A virus hemagglutinin subtype (H16) obtained from black-headed gulls[END_REF]) is related to the distances between the C-atoms, in the reference and target peptide planes, and the local axis of screw motion. The radius ρ is a measure of the local curling of the backbone conformation. In a flat backbone conformation, like an extended β-strand, ρ is close to zero ; when the local backbone conformation is curled, as in the case of α-helices and β-turns, ρ increases. The maximum observed values usually do not exceed 0.3[START_REF] Kneller | Efficient characterization of protein secondary structure in terms of screw motions[END_REF][START_REF] Calligari | ScrewFit: a novel approach for continuum protein secondary structure assessments[END_REF].

  ), give properties intrinsically related to the pair of consecutive peptide planes and should not be associated only to one residue. This fact must be remembered when interpreting the ScrewFit profiles shown in this work.The comparison of ScrewFit profiles of different NAs is performed through the punctual differences along the sequence axis. This type of comparison is possible only because of the absence of insertions or deletions in the multi-alignment of the corresponding sequences. In the opposite case, comparison is possible only on limited regions where this condition is still verified. The positive differences in the radius of screw motion (ρ) and the orientational distance ( Δ Ω ), indicate a change toward a more curled conformation of the backbone, whereas any deviation from zero of the straightness (σ) is caused by a kink (if located in one or two residues) or a smoother curvature (if more than two residues are involved) in one of the two main chains.

  appears to interact with oseltamivir in group-1 and group-2 NAs in a similar fashion, its geometry is slightly affected by the open-to-close conformational transition[START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF]. The variation of the ScrewFit parameters (see Fig.2) suggests that the NA1 open-to-closed transition parallels with the shift of the Glu119 residue towards Arg118; this could induce a small bending of the Phe115-Cys124 β-strand as indicated by the value of straightness parameter between Ile117 and Glu119.In region A of NA1, the oseltamivir-induced structural difference(s) distinguishing the open to the closed conformations entail(s) the reduction of the helical radius of several residues and the increase of the orientational distance (see Fig.2). In the closed conformation, the peptide planes of residues Leu140-Arg152 tend to be more distant from the local rotation axis of the screw motion, indicating a change from a flat to a curled conformation. This result can be inferred from the reduction of the radius of the screw motion (ρ) and the corresponding change in the orientational distance ( Δ Ω ). The two parameters show that consecutive plans attain different orientations, either lowering their orientational distance (residues Arg130-Ser145) or increasing the latter (region Gly147-Ser160), which demand for a global rearrangement around the common axis of screw motion. Such a conformational change is also reflected by a remarkable change of the local rotation axis in the sub-region Leu140-Ser145, as shown by the value of the straightness difference which is grater than one, this indicating that the transition between open and closed conformation induces a change of the straightness from a non-negative value to a non-positive one (see Fig.2). The latter reads as a change of about π/2 in the direction of the local rotation axis. The differences of the ScrewFit parameters match with the moving of the subregion Ser145-Arg152 to the center of the active site[START_REF] Russell | The structure of H5N1 avian influenza neuraminidase suggests new opportunities for drug design[END_REF] (Fig.3B).Changes of theScrewFit parameters of region B are similar to those observed in region A, the backbone attaining a more packed conformation around Asn247 as revealed by A C C E P T E D M A N U S C R I P T the positive change in the straightness (close to one), which indicates a kink into the local axis of rotation between residues Asn247 and Gly248. Nevertheless, it is worth to note that here the screw motion helix radius decreases from the open to the closed conformation meaning that the peptide planes assume similar orientations. The extent of this change (~ 0.1) shows that peptide planes between residue Asn247 and Gly248 adopt a β-turn-like configuration in the closed oseltamivir-bound NA1. The variation in the NA1 structure involves Ala250,

  ScrewFit plot for the NA1 closed conformation minus ligand-free NA9. Although very similar, our results show that in regions A, B, and C both the open and the closed configuration of NA1 are characterized by variations of the ScrewFit parameters with respect to the NA9 configuration.

Fig. 1 .

 1 Fig. 1. A tri-peptide with two peptide bonds in the extended conformation, where the symbol

Fig. 2 .

 2 Fig. 2. Screw-Fit parameter differences between the free and oseltamivir-bound open NA1

Fig. 3 .

 3 Fig. 3. Comparison between the free, open and closed structures of NA1. ScrewFit

Fig. 4 .Fig. 5 .

 45 Fig. 4. ScrewFit parameter differences of the ligand-free minus the oseltamivir-bound NA9

Figure 1 A

 1 Figure 1

Figure 3 A

 3 Figure 3

Acknowledgments

Authors wish to thank Dr. Manuela Emili and Dr. Angelo Merante for graphical assistance.

Figure 5

(cyan). The oseltamivir is reported in ball-and-stick and colored according to atom type.

Picture was drawn by means of UCSF Chimera [START_REF] Pettersen | UCSF Chimera -a visualization system for exploratory research and analysis[END_REF].