Andrea Braides
email: braides@mat.uniroma2.it

Antoine Gloria
email: gloria@cermics.enpc.fr

EXACT BOUNDS ON THE EFFECTIVE BEHAVIOUR OF A CONDUCTING DISCRETE POLYCRYSTAL

Keywords: polycrystals, discrete energies, effective properties, composite materials AMS subject classifications. 77Q15, 74S05, 35J60

In a recent paper by Braides and Francfort, the problem of the characterization of the overall properties of lattice energies describing networks with arbitrary mixtures of two types of linear conductors has been addressed in a two-dimensional setting. In this paper we investigate the connection between that discrete optimization process and the theory of bounds for mixtures of continuum energies, for which the choice of the relationships between the different phases of the mixture is unusual and leads to remarkably simple results in terms of G-closure.

1. Introduction. In a paper by Braides and Francfort [START_REF] Braides | Bounds on the effective behaviour of a square conducting lattice[END_REF], the problem of the characterization of the overall properties of lattice energies describing networks with arbitrary mixtures of two types of linear conductors has been addressed. In a twodimensional setting, with the notation of Γ-convergence, that problem translates into the limit analysis as h → 0 of discrete energies of the type

E h (u) = 1 2 (i,j)∈N h c h ij (u i -u j ) 2 u : N h → R,
where N h denotes the set of nearest neighbours (i.e., pairs (i, j) such that |i -j| = h) on a portion of a square lattice hZ 2 ∩ Ω, Ω being a regular open subset of R 2 , and

c h ij ∈ {α, β}
are chosen arbitrarily. It should be noted that in the extreme cases c h ij ≡ α and c h ij ≡ β the functionals above are simply discretizations of the Dirichlet integral with the corresponding coefficient. Since the analogue of the 'localization principle' holds in this discrete setting, the Γ-limit problem translates in a sort of 'discrete G-closure problem', where bounds on a matrix given by a 'discrete homogenization formula' must be exhibited. The results in [START_REF] Braides | Bounds on the effective behaviour of a square conducting lattice[END_REF] show that all possible limits of E h have the form

F (u) = Ω A(x)∇u, ∇u dx u ∈ H 1 (Ω),
with A(x) belonging to a set of matrices determined only by the local limit proportion θ(x) at x of 'α-connections'. Such sets can be compared with those obtained by limits of mixtures of isotropic energies; i.e., of functionals of the form with Ω ǫ arbitrary measurable subsets of Ω (χ Ωǫ denotes the corresponding characteristic function). These limits have the same form F with A(x) now belonging to a set whose form depends on the limit volume fraction θ(x) of the α-phase; i.e., on the value at x of the weak * -limit of χ Ωǫ (see Murat and Tartar [START_REF] Murat | H-convergence[END_REF][START_REF] Murat | Séminaire d'Analyse Fonctionnelle et Numérique[END_REF] and Cherkaev and Lurie [START_REF] Lurie | Exact estimates of conductivity of mixtures composed of two materials taken in prescribed proportion (plane problem)[END_REF]). In Figure 1.1 the two problems are compared in terms of the possible eigenvalues of A(x) for some θ(x) < 1/2: the shaded sets represent all possible pairs of eigenvalues in the two cases, the smaller set being the 'continuous' one (see Tartar [START_REF] Tartar | Estimation de coefficients homogénéisés[END_REF][START_REF] Tartar | Estimations fines de coefficients homogénéisés[END_REF]). This comparison shows that the discrete energies cannot be simply interpreted as a mixture of the two 'extremal' isotropic energies.

E ǫ (u) = Ω αχ Ωǫ + β(1 -χ Ωǫ ) |∇u| 2 dx u ∈ H 1 (Ω),
The discrete functionals can be heuristically interpreted as 'anisotropic' continuous ones in the following way, pictured in Figure 1.2. We introduce a triangulation with underlying lattice hZ 2 , and to every triangle T with a horizontal side corresponding to a connection with value c h ij = c h and a vertical side corresponding to a connection with c h ij = c v we associate the matrix

A h (x) = c h 0 0 c v for x ∈ T.
In Figure 1.2, on the left, 'thick connections' correspond to the coefficient β (the other ones to α), while on the right on black triangles we take A h (x) = βI, on white triangles A h (x) = αI, and on the remaining triangles the anisotropic conductivities.

Note that A h takes only four possible values. Chosing ũh ∈ H 1 (Ω) as the piecewise affine interpolation of u h defined on N h , we have

E h (u h ) = E h (ũ h ) + o(1),
(the last term deriving from an asymptotically negligible boundary term; see e.g. [START_REF] Alicandro | Representation result for continuum limits of discrete energies with superlinear growth[END_REF]) where

E h (u) = Ω A h (x)∇u, ∇u dx u ∈ H 1 (Ω).
Note that the infimum of E h on H 1 (Ω) is in general strictly lower than the corresponding one for E h provided compatible boundary conditions are given. On the other hand, any continuous configuration can be approximated by a discrete configuration (let simply think of the discretization of the Dirichlet integral). Therefore, the Γ-closure of the discrete energies E h contains the Γ-closure of the continuous energies E h . Whether this inclusion is strict or not is unclear. However, we will show that for polycristalline structures the two closures coincide. Continuous energies as E h can also be thought of as being a polycrystalline mixture of three conducting materials by regarding one anisotropic value of A h as the rotation of the other one, so that we may rewrite

E h (u) = Ω A h (x)R h (x)∇u(x), R h (x)∇u(x) dx, u ∈ H 1 (Ω), with A h (x) taking one of the values α 0 0 α , β 0 0 β , or α 0 0 β , (1.1)
and the rotation matrix R h (x) being either the identity or the rotation by ninety degrees. A more natural setting for the continuous energy is obtained by removing the geometric restrictions due to the lattice symmetries, and determine all the possible conductivity tensors that can be obtained as possible Γ-limits of energies

F ǫ (u) = Ω A ǫ (x)R ǫ (x)∇u(x), R ǫ (x)∇u(x) dx, (1.2)
with A ǫ and R ǫ families of matrices, with A ǫ taking the three values above, and R ǫ arbitrary rotations. If one adds the restriction that R ǫ be piecewise-constant, then we may interpret this continuous energy as a limit of discrete polycrystals, that is a polycrystal obtained by a discrete to continuum process. This discrete polycrystalline material is composed of grains of diameter much larger than the lattice spacing h, at whose interior we consider a network structure whose lattice orientation is described by the rotation matrix R ǫ (x), as depicted in Figure 1.3 (the black and grey subdomains represent isotropic parts of the energy whereas the other subdomains correspond to underlying lattices oriented according to the arrows). In this setting, the characterization of the Γ-closure of such discrete polycrystals is complete since the additional invariance by rotation of the Γ-limit allows us to focus on diagonal matrices, which have been fully described in [START_REF] Braides | Homogenization of Multiple Integrals[END_REF].

The scope of the present paper is twofold: to prove that the bounds obtained by considering piecewise-constants rotations R ǫ still hold for general rotations, and that these bounds -which are optimal in the discrete setting -are also attained in this continuous setting. This amounts to proving that the Γ-closures of discrete and continuous polycrystals are the same. In what follows, these two types of polycrystal will be refered to as 'discrete' since the Γ-closures are related to the very special relationship between the phases, that is inherited from the discrete description.

From a modeling viewpoint, these energies can equivalently be thought of as describing a bilayered polycrystal made of two types of layers (namely A and B). Each layer has a conductivity tensor Cν ⊗ ν, where ν ∈ S 1 (R 2 ), C = α for type A and C = β for type B, with β > α > 0. Each simple crystal is a bilayer C 1 ν 1 ⊗ν 1 +C 2 ν 2 ⊗ν 2 with C i ∈ {α, β} obtained by the superposition of two simple layers in orthogonal directions ν 1 ⊥ ν 2 , such that the conductivity tensor of the crystal is of the form above (1.1). The association of rotated simple crystals gives rise to the polycrystal under investigation.

It is convenient to rewrite the energies F ǫ in (1.2) in a different form, as

F ǫ (u) = α E ǫ 1 |∇u| 2 dx + β Ω\E ǫ 2 |∇u| 2 dx + E ǫ 2 \E ǫ 1 α 0 0 β R ǫ (x)∇u(x), R ǫ (x)∇u(x) dx, (1.3)
where

E ǫ 1 ⊂ E ǫ 2 .
The set E ǫ 1 is the part of the polycrystal with isotropic conductivity α, Ω \ E ǫ 2 is the part of the polycrystal with isotropic conductivity β, and E ǫ 2 \ E ǫ 1 is where each single crystal has an anistropic conductivity of type α 0 0 β , up to a rotation. In order to describe the limit of F ǫ , up to the extraction of a subsequence, we may

suppose that χ E ǫ 1 ⇀ θ 1 and χ E ǫ 2 ⇀ θ 2 weakly * in L ∞ (Ω)
, and we define θ =

θ 1 + θ 2 2 .
In the notation above the value θ(x) represents the fixed limit local proportion of layers of type A at the point x, and is the analog of the local proportion of α-connections in the discrete setting. With the classical notation of G-closure, we would denote by θ α the volume fraction of isotropic conductivity α, by θ β the volume fraction of isotropic conductivity β and by θ α,β = 1θ αθ β the volume fraction of anisotropic conductivity. The problem we are interested in is not the classical G-closure of a three phase mixture, which is a very tough issue. Actually our discrete setting imposes an additional relationship between θ α , θ α and θ α,β , namely θ α + θ α,β 2 = θ, or equivalently

θ β + θ α,β 2 = 1 -θ,
where θ is fixed by the local proportion of α-connections. Our aim is then to characterize the Γ-limits for all values of θ ∈ [0, 1]. Surprisingly, this additional relationship trivializes the associated G-closure problem, which can be interpreted as a 'constrained' three-phase mixture G-closure problem.

The paper is organized as follows. In the second section, we give a rigorous formulation of the result that we prove in two steps: first we derive bounds on the eigenvalues of A 0 (x), then we show these bounds to be optimal. We perform a derivation of these bounds following a purely variational approach to highlight how our parameter θ intervenes in their computations, so that the optimality of those bounds is simply obtained by making the previous reasoning sharp. This proof is direct and self-contained, and can be compared with the one in the discrete setting in [START_REF] Braides | Bounds on the effective behaviour of a square conducting lattice[END_REF]. In the third section, we generalize the result to the uniaxial three-dimensional case. The bounds we get in the continuous case are indeed the 'trivial bounds' usually obtained in the continuum theory of bounds (for a general exposition of results connected to the computation of bounds for the effective properties of composites we refer to the monograph by Milton [START_REF] Milton | The Theory of Composites[END_REF]), which may be surprising at the first glance since we deal with a three-phase mixture. The last section is dedicated to the reformulation and interpretation of the present results in terms of a constrained G-closure problem, where we recall the more classical proof of the trivial bounds and generalize the results.

Main result and derivation of bounds in two dimensions.

2.1. Optimal bounds. Having set the notation in the Introduction, we may describe the Γ-limits of functionals F ǫ as follows.

Theorem 1. Let F ǫ be defined by

(1.3), let χ E ǫ 1 ⇀ θ 1 and χ E ǫ 2 ⇀ θ 2 weakly * in L ∞ (Ω), and let θ = θ 1 + θ 2 2 .
Up to the extraction of a (further) subsequence F ǫ Γ-converges to an energy of the form

F 0 (u) = Ω A 0 (x)∇u, ∇u dx (2.1)
on H 1 (Ω), where for almost all x ∈ Ω we have A 0 (x) ∈ H 2 d (θ(x)), and for fixed θ ∈ [0, 1] the set H 2 d (θ) is defined as the set of symmetric 2 × 2 matrices whose eigenvalues λ 1 , λ 2 satisfy the bounds

α ≤ λ 1 , λ 2 ≤ β, λ 1 + λ 2 ≤ 2 (θα + (1 -θ)β) , 1 λ 1 + 1 λ 2 ≤ 2 θ α + 1 -θ β (2.2) Conversely, for all measurable θ and A 0 such that A 0 (x) ∈ H 2 d (θ(x)), there exist E ǫ i as above such that F ǫ Γ-converges to F 0 in (2.

1).

Theorem 1 states that the effective overall conductivity of the polycrystal can be characterized in terms of the local proportion of layers of type A. As motivated in the Introduction, these bounds correspond to the bounds obtained for a square conducting lattice in [START_REF] Braides | Bounds on the effective behaviour of a square conducting lattice[END_REF] for diagonal matrices, or for a discrete polycrystal. With this observation in mind we have added the subscript 'd' for discrete in the notation above. In terms of those bounds it must be noted that the optimal micro-geometries correspond to lattices locally oriented in the directions of the eigenvectors of the matrix A 0 (x).

For fixed θ ∈ [0, 1] we will denote by H 2 hom (θ) the set of all homogenized matrices obtained by the homogenization of a periodic polycrystal as above with underlying volume fraction θ of A-layers; i.e., of all matrices B that satisfy the equality

Bξ, ξ = inf α E1 |ξ + ∇u| 2 dy + β (0,1) 2 \E2 |ξ + ∇u| 2 dy + E2\E1 α 0 0 β R(y)(ξ + ∇u(y)), R(y)(ξ + ∇u(y)) dy : u ∈ H 1 # ((0, 1) 2 ) , (2.3
) for some rotation R(y), and subsets

E 1 ⊂ E 2 ⊂ (0, 1) 2 , with |E 1 | + |E 2 | = 2θ.
Here and after H 1 # ((0, 1) 2 ) denotes the space of H 1 loc (R 2 ) functions that are 1-periodic in the coordinate directions. Note that the energy Ω B∇u, ∇u dx corresponding to such B is a particular Γ-limit of a family of the type (1.3), with E ǫ i = ǫE i (E i are extended by 1-periodicity) and R ǫ (y) = R(y/ǫ). The triplet (E 1 , E 2 , R) is called the (micro-)geometry of the corresponding periodic polycrystal.

First, note that by a general compactness result for the Γ-convergence of quadratic forms ( [START_REF] Braides | Homogenization of Multiple Integrals[END_REF], [START_REF] Braides | Γ-convergence for Beginners[END_REF], [START_REF] Maso | An Introduction to Γ-convergence[END_REF], [START_REF] Braides | A handbook of Γ-convergence[END_REF] Theorem 4.13) we obtain that, up to subsequences, F ǫ Γconverges to some F 0 represented as in (2.1) for some A 0 (x). By the well known localization principle (see [START_REF] Braides | A handbook of Γ-convergence[END_REF] Section 5.4.1 for a simple version, or [START_REF] Raitums | On the local representation of G-closure[END_REF] for a slightly more general one), we have that for almost every x ∈ Ω and for all ξ ∈ R 2 ,

A 0 (x) ∈ H 2 hom (θ(x) + o(1)) + o(1); (2.4)
i.e., for every η > 0 there exist R, E 1 , and

E 2 , with ||E 1 | + |E 2 | -2θ(x)| < η, such that if B is as above then |B -A 0 (x)| ≤ o(1) as η → 0.
In the rest of the section we will prove that H 2 hom (θ) = H 2 d (θ). From this inequality, the continuity of the bounds in (2.2), and the closure of

H 2 d (θ), it then follows that A 0 (x) ∈ H 2 d (θ(x)) = H 2 hom (θ(x)) for a.a.
x ∈ Ω and the first part of Theorem 1. We will therefore fix E 1 , E 2 , and R as in (2.3), and introduce a rotation angle φ. so that we can write R as R(y) = cos(φ(y)) sin(φ(y)) sin(φ(y)) cos(φ(y)) .

Note that the set H 2 hom (θ) is invariant by multiplication by rotations. Since it is also composed of symmetric matrices, this set is consequently characterized by the set of all possible eigenvalues of such matrices. The bounds are therefore derived in the plane of the eigenvalues (λ 1 , λ 2 ) of matrices in H 2 hom (θ). 2.2. Simple bounds. Proposition 2.1. For i = 1, 2 we have α ≤ λ i ≤ β. These easy bounds follow immediately by the pointwise estimate α Id ≤ A ǫ (y) ≤ β Id.

Bounds from arithmetic means.

Proposition 2.2. We have

λ 1 + λ 2 ≤ 2 (αθ + (1 -θ)β) .
(2.5)

Proof. Take ξ = e i and the test function u = 0 in (2.3). We then have for i = 1, 2

λ i ≤ α|E 1 | + β(1 -|E 2 |) + E2\E1 α 0 0 β R(y)e i , R(y)e i dy (2.6)
In addition, since {e 1 , e 2 } is an orthonormal basis of R 2 and R(y) is a rotation, for almost every y ∈ E 2 \ E 1 we have

α 0 0 β R(y)e 1 , R(y)e 1 + α 0 0 β R(y)e 2 , R(y)e 2 = α + β. (2.7)
Indeed the first term of (2.7) reads α cos 2 φ + β sin 2 φ and the second term β cos 2 φ + α sin 2 φ. The combination of (2.6) and (2.7) gives

λ 1 + λ 2 ≤ 2α|E 1 | + 2β(1 -|E 2 |) + (|E 2 | -|E 1 |)(α + β), (2.8)
which yields (2.5)

using θ = |E 2 | + |E 1 | 2 .
2.4. Bounds from harmonic means. Proposition 2.3. We have

1 λ 1 + 1 λ 2 ≤ 2 θ α + 1 -θ β .
(2.9)

The rest of the section is devoted to proving this bound, relying on two arguments. The first one is a slicing argument, which is based on the simple inequality |∇u| 2 ≥ (∂ i u) 2 and allows to reduce to one-dimensional problems. The second tool is some convexity inequality.

Slicing argument. We first prove the following inequalities:

λ 1 ≥ 1 0 1 α |E 1 1 (y 2 )| + 1 β (1 -|E 1 2 (y 2 )|) + E 1 2 (y2)\E 1 1 (y2) α sin 2 φ(y) + β cos 2 φ(y) αβ dy 1 -1 dy 2 , λ 2 ≥ 1 0 1 α |E 2 1 (y 1 )| + 1 β (1 -|E 2 2 (y 1 )|) + E 2 2 (y1)\E 2 1 (y1) β sin 2 φ(y) + α cos 2 φ(y) αβ dy 2 -1 dy 1 , (2.10)
where

E 1 1 (y 2 ) = {y 1 ∈ (0, 1), χ E1 (y 1 , y 2 ) = 1}, E 1 2 (y 2 ) = {y 1 ∈ (0, 1), χ E2 (y 1 , y 2 ) = 1}, E 2 1 (y 1 ) = {y 2 ∈ (0, 1), χ E1 (y 1 , y 2 ) = 1}, E 2 2 (y 1 ) = {y 2 ∈ (0, 1), χ E2 (y 1 , y 2 ) = 1},
and χ E denotes the characteristic function of the set E. These sets satisfy

1 0 |E 1 1 (y 2 )|dy 2 = 1 0 |E 2 1 (y 1 )|dy 1 = θ 1 and 1 0 |E 1 2 (y 2 )|dy 2 = 1 0 |E 2 2 (y 1 )|dy 1 = θ 2 .
To get (2.10), we first prove the following result.

Proposition 2.4. For all (u 1 , u 2 ) ∈ R 2 and φ ∈ (0, 2π), we have

α 0 0 β R(φ)(u 1 e 1 + u 2 e 2 ), R(φ)(u 1 e 1 + u 2 e 2 ) ≥ u 2 1 αβ α sin 2 φ + β cos φ , (2.11)
α 0 0 β R(φ)(u 1 e 1 + u 2 e 2 ), R(φ)(u 1 e 1 + u 2 e 2 ) ≥ u 2 2 αβ β sin 2 φ + α cos φ , (2.12
) where {e 1 , e 2 } is the canonical basis of R 2 . Proof. It is enough to prove (2.11). Denote by χ φ (u 1 , u 2 ) the left-hand side of (2.11). A direct calculation shows that

χ φ (u 1 , u 2 ) = α(u 1 cos φ -u 2 sin φ) 2 + β(u 1 sin φ + u 2 cos φ) 2 . At fixed u 1 , v → χ φ (u 1 , v
) is quadratic and convex. This function attains its minimum on R, at

v 2 = u 1 α -β α sin 2 φ + β cos 2 φ
sin φ cos φ.

We then have

χ φ (u 1 , v 2 ) = u 2 1 α cos 2 φ + β sin 2 φ - (β -α) 2 sin 2 φ cos 2 φ α sin 2 φ + β cos 2 φ = u 2 1 αβ cos 2 φ + sin 2 φ 2 α sin 2 φ + β cos 2 φ = u 2 1 αβ α sin 2 φ + β cos 2 φ as desired.
We now integrate in the e 2 direction. If u is such that ∇u = e 1 + ∇w, w ∈ H 1 # ((0, 1)), we then have

λ 1 = 1 0 α E 1 1 (y2) |∇u| 2 dy 1 + β (0,1)\E 1 2 (y2) |∇u| 2 dy 1 + E 1 2 (y2)\E 1 1 (y2) α 0 0 β R(y)∇u, R(y)∇u dy 1 dy 2 ≥ 1 0 α E 1 1 (y2) ∂ 1 u 2 dy 1 + β (0,1)\E 1 2 (y2) ∂ 1 u 2 dy 1 + E 1 2 (y2)\E 1 1 (y2) dy 1 ∂ 1 u 2 αβ α sin 2 φ(y) + β cos 2 φ(y) dy 2 ≥ 1 0 inf α E 1 1 (y2) v ′ (y 1 ) 2 dy 1 + β (0,1)\E 1 2 (y2) v ′ (y 1 ) 2 dy 1 + E 1 2 (y2)\E 1 1 (y2) v ′ (y 1 ) 2 αβ α sin 2 φ(y) + β cos 2 φ(y) dy 1 : v ∈ H 1 0 ((0, 1)) dy .
(2.13) The infimum in the integral of (2.13) can be explicitly computed by using the following well-known result.

Proposition 2.5. Let f ∈ L ∞ ((0, 1)) be such that inf f > 0, then

inf 1 0 f (y)v ′ (y) 2 : v ∈ H 1 (0, 1), v(0) = 0, v(1) = 1 = 1 0 dy f (y) -1 (2.14)
Inserting formula (2.14) in (2.13), we obtain

λ 1 ≥ 1 0 1 α |E 1 1 (y 2 )| + 1 β (1 -|E 1 2 (y 2 )|) + E 1 2 (y2)\E 1 1 (y2) α sin 2 φ(y) + β cos 2 φ(y) αβ dy 1 -1 dy 2 .
(2.15)

The same proof holds for λ 2 .

To conclude the proof of the bounds, we use some convexity properties.

A convexity inequality. Proposition 2.6. Let b > a > 0, c > 0, γ > 0 and

f ∈ L 1 ((a, b), R + ), then b a dy c + f (y) ≥ (b -a) 2 b a (c + f (y))dy (2.16)
This result just amounts to stating that the harmonic mean is lower than the arithmetic mean, and it is immediately obtained by Hölder's inequality. Applying Proposition 2.6 to (2.10), for a = 0, b = 1 and c = 1/α, we obtain

λ 1 ≥ 1 0 1 α |E 1 1 (y 2 )| + 1 β (1 -|E 1 2 (y 2 )|) + E 1 2 (y2)\E 1 1 (y2) α sin 2 φ(y) + β cos 2 φ(y) αβ dy 1 dy 2 -1 = θ 1 α + 1 -θ 2 β + E2\E1 α sin 2 φ(y) + β cos 2 φ(y) αβ dy -1 λ 2 ≥ 1 0 1 α |E 2 1 (y 1 )| + 1 β (1 -|E 2 2 (y 1 )|) + E 2 2 (y1)\E 2 1 (y1)
β sin 2 φ(y) + α cos 2 φ(y) αβ dy 2 dy 1

-1 = θ 1 α + 1 -θ 2 β + E2\E1 β sin 2 φ(y) + α cos 2 φ(y) αβ dy -1
(2.17)

Therefore we have

1 λ 1 + 1 λ 2 ≤ 2 θ 1 α + 2 1 -θ 2 β + E2\E1 α cos 2 φ + β sin 2 φ αβ + β cos 2 φ + α sin 2 φ αβ dy = 2 θ 1 α + 2 1 -θ 2 β + (θ 2 -θ 1 ) 1 α + 1 β
which exactly reads as (2.9).

3. Optimality of the bounds. We now prove the last statement of Theorem 1. By approximation with piecewise-constant functions, it is sufficient to consider the case of θ and A 0 constant; i.e., it will be sufficient to show that for all symmetric A 0 with eigenvalues satisfying (2.2) we find a suitable geometry such that each A 0 ξ, ξ is represented as a minimum problem as in (2.3). By rotational invariance, it suffices to consider diagonal A 0 . Furthermore, it is not restrictive to treat the case of eigenvalues belonging to the boundary of the set defined by (2.3) only (the 'extremal' cases). In fact, using a well-known "lamination formula" (see [START_REF] Tartar | Estimations fines de coefficients homogénéisés[END_REF][START_REF] Tartar | Estimation de coefficients homogénéisés[END_REF] Proposition 3) we may easily construct geometries for a generic diagonal A 0 using those 'extremal geometries'.

We will consider polycrystals (0, 1) 2 defined by A(y) = a 1 (y) 0 0 a 2 (y) for y ∈ (0, 1) 2 , and the corresponding homogenized A 0 , defined by

A 0 ξ, ξ = inf (0,1) 2 A(y)(ξ + ∇u(y)), ξ + ∇u(y) , u ∈ H 1 # ((0, 1) 2 ) (3.1)
for all ξ ∈ R 2 . Note that in terms of the micro-geometry of the polycrystal the rotation R is either the identity or a rotation of ninety degrees.

3.1. Arithmetic bound. Let θ 1 ∈ (0, 1) and θ 2 ∈ (0, 1). We define A(y) by We will then prove that

a 1 (y) = α if y 2 ≤ θ 1 , a 1 (y) = β if y 2 > θ 1 , and a 2 (y) = α if y 1 ≤ θ 2 , a 2 (y) = β if y 1 > θ 2 . α 0 0 α 0 β β 0 0 α β 0 0 β α 0 ( ) ( ) ( ) ( )
A 0 = ā(θ 1 ) 0 0 ā(θ 2 )
, where ā(λ) = αλ + β(1λ).

Testing with u = 0 in (3.1), we obtain A 0 ξ, ξ ≤ ā(θ 1 )ξ 2 1 + ā(θ 2 )ξ 2 2 . To prove the converse inequality we use the slicing method, obtaining

(0,1) 2 A(y)(ξ + ∇u(y)), ξ + ∇u(y) = 1 0 θ1 0 α(ξ 1 + ∂ 1 u(y 1 , y 2 )) 2 dy 2 dy 1 + 1 0 1 θ1 β(ξ 1 + ∂ 1 u(y 1 , y 2 )) 2 dy 2 dy 1 + 1 0 θ2 0 α(ξ 2 + ∂ 2 u(y 1 , y 2 )) 2 dy 1 dy 2 + 1 0 1 θ2 β(ξ 2 + ∂ 2 u(y 1 , y 2 )) 2 dy 1 dy 2 . (3.2)
Switching the order of the integrals and minimizing each term of the right-hand side of (3.2) separately yield the desired inequality, and thus the equality

A 0 ξ, ξ = ā(θ 1 )ξ 2 1 + ā(θ 2 )ξ 2 2 .
If θ 2 + θ 1 = 2θ then ā(θ 1 ) + ā(θ 2 ) = αθ + β(1θ), which is the arithmetic bound.

Harmonic bound.

For the harmonic bound, let us consider a 1 (y) = α if

y 1 ≤ θ 1 , a 1 (y) = β if y 1 > θ 1 , and a 2 (y) = α if y 2 ≤ θ 2 , a 2 (y) = β if y 2 > θ 2 . α 0 0 α 0 β ( ) ( ) ( ) ( ) α 0 0 β β 0 β 0 0 α θ θ 2 1 0 1 1 0 Fig. 3.2.
Optimal geometry for the harmonic bound

We will then prove that

A 0 = a(θ 1 ) 0 0 a(θ 2 )
, where a(λ

) = λ α + 1 -λ β -1
.

Let u 1 : (0, 1) → R be the one-dimensional minimizer of

inf θ1 0 α(ξ 1 + v ′ (z)) 2 dz + 1 θ1 β(ξ 1 + v ′ (z)) 2 dz : v(0) = v(1) = 0 = a(θ 1 )ξ 2 1 ,
and u 2 : (0, 1) → R be the one-dimensional minimizer of

inf θ2 0 α(ξ 2 + v ′ (z)) 2 dz + 1 θ2 β(ξ 2 + v ′ (z)) 2 dz : v(0) = v(1) = 0 = a(θ 2 )ξ 2 2 .
Let u(y) = u 1 (y 1 ) + u 2 (y 2 ) be a test function in (3.1). We then have

A 0 ξ, ξ ≤ a(θ 1 )ξ 2 1 + a(θ 2 )ξ 2 2 .
Once more we prove the converse inequality using the slicing argument.

(0,1) 2

A(y)(ξ + ∇u(y)), ξ + ∇u(y) = 1 0 θ1 0 α(ξ 1 + ∂ 1 u(y 1 , y 2 )) 2 dy 1 + 1 θ1 β(ξ 1 + ∂ 1 u(y 1 , y 2 )) 2 dy 1 dy 2 + 1 0 θ2 0 α(ξ 2 + ∂ 2 u(y 1 , y 2 )) 2 dy 2 + 1 θ2 β(ξ 2 + ∂ 2 u(y 1 , y 2 )) 2 dy 2 dy 1 (3.
3) We obtain the desired inequality by optimizing pointwise the integrands of the right handside of (3.3).

If θ 2 + θ 1 = 2θ then 1 a(θ 1 ) + 1 a(θ 2 ) = θ α + 1 -θ β
, which is the harmonic bound.

3.3. Optimal bounds. The simple bounds λ i = α, β can always be reached by layering. Depending on θ, the shape of the optimal bounds may vary (see Figure 5 in [START_REF] Braides | Bounds on the effective behaviour of a square conducting lattice[END_REF], and the details therein). Once the arithmetic and harmonic bounds are shown to be attained, the convex set delimited by these bounds and the simple bounds can also be attained by layering.

Extension to higher dimension

. Theorem 1 has a natural counterpart in higher dimension. We only consider the physically meaningful three-dimensional case, and especially the uniaxial basic crystal. In this case we deal with the integral functionals

F ǫ (u) = α E ǫ 1 |∇u| 2 + β Ω\E ǫ 2 |∇u| 2 + E ǫ 3   α 0 0 0 β 0 0 0 β   R ǫ (x)∇u(x), R ǫ (x)∇u(x) dx + (E ǫ 2 \E ǫ 1 )\E ǫ 3   β 0 0 0 α 0 0 0 α   R ǫ (x)∇u(x), R ǫ (x)∇u(x) dx, (4.1) where E ǫ 1 ⊂ E ǫ 2 , E ǫ 3 ⊂ E ǫ 2 \ E ǫ
1 and R ǫ is a rotation. As in the two-dimensional case, in order to describe the Γ-limit of F ǫ , up to the extraction of a subsequence, we may suppose that χ

E ǫ 1 ⇀ θ 1 , χ E ǫ 2 ⇀ θ 2 , and χ E ǫ 3 ⇀ θ 3 weakly * in L ∞ (Ω), and we define θ = θ 1 + 2θ 2 -θ 3 3
, again the value θ(x) representing the fixed limit local proportion of layers of type A at the point x.

Theorem 2. Let F ǫ be defined by (4.1), let χ

E ǫ i ⇀ θ i weakly * in L ∞ (Ω), and let θ = θ 1 + 2θ 2 -θ 3 3 .
Up to the extraction of a (further) subsequence F ǫ Γ-converges to an energy of the form

F 0 (u) = Ω A 0 (x)∇u, ∇u dx (4.2)
on H 1 (Ω), where for almost all x ∈ Ω we have A 0 (x) ∈ H 3 d (θ(x)), where for fixed θ ∈ [0, 1] the set H 3 d (θ) is defined as the set of symmetric 3 × 3 matrices whose eigenvalues λ 1 , λ 2 , λ 3 satisfy the bounds

α ≤ λ 1 , λ 2 , λ 3 ≤ β, λ 1 + λ 2 + λ 3 ≤ 3ā(θ), 1 λ 1 + 1 λ 2 + 1 λ 3 ≤ 3 a(θ) (4.3)
Conversely, for all measurable θ and A 0 such that A 0 (x) ∈ H 3 d (θ(x)), there exist E ǫ i as above such that F ǫ Γ-converges to F 0 in (4.2). We now adapt Sections 2 and 3 to prove Theorem 2.

Since the anisotropic conductivity matrices only have two eigenvalues α and β, in order to diagonalize the matrix it suffices to determine the eigenvector associated to the eigenvalue of multiplicity 1, the orthogonal subspace being an eigensubspace associated to the other eigenvalue. Therefore, we only need two angles for the rotations which can be written as

R(φ, ψ) =   sin ψ cos φ cos ψ sin φ cos ψ cos ψ -cos φ sin ψ -sin φ sin ψ 0 -sin φ cos φ   .
The calculations are more involved in the three-dimensional case. We leave the details to the interested reader and only give the intermediate results.

Bounds from arithmetic means. Direct calculations give

  α 0 0 0 β 0 0 0 β   R(φ, ψ)e 1 , R(φ, ψ)e 1 = α sin 2 ψ + β cos 2 ψ   α 0 0 0 β 0 0 0 β   R(φ, ψ)e 2 , R(φ, ψ)e 2 = α cos 2 φ cos 2 ψ + β(cos 2 φ sin 2 ψ + sin 2 φ)   α 0 0 0 β 0 0 0 β   R(φ, ψ)e 3 , R(φ, ψ)e 3 = α sin 2 φ cos 2 ψ + β(sin 2 φ sin 2 ψ + cos 2 φ),
whose sum is α + 2β. The calculation for the other anisotropic conductivity matrix yields β + 2α. Therefore,

λ 1 + λ 2 + λ 3 ≤ 3θ 1 α + 3(1 -θ 2 )β + θ 3 (α + 2β) + (θ 2 -θ 1 -θ 3 )(2α + β) = 3(αθ + (1 -θ)β),
which is the desired arithmetic bound.

4.2. Bounds from harmonic means. For convenience, we will make use of c for cos, s for sin, and v = (v 1 , v 2 , v 3 ) in the sequel. We set

H(v 1 , v 2 , v 3 ) :=   α 0 0 0 β 0 0 0 β   R(φ, ψ)v, R(φ, ψ)v = α(v 1 s ψ + v 2 c φ c ψ + v 3 s φ c φ ) 2 + β(v 1 c ψ -v 2 c φ s ψ -v 3 s φ s ψ ) 2 +β(-v 2 s φ + v 3 c φ ) 2
As in the two-dimensional case, we optimize H on R 2 , one variable being fixed. These minimizers exist since H is a symmetric positive definite quadratic form. The first derivatives of H with respect to v 1 , v 2 , v 3 read as

∂ 1 H(v 1 , v 2 , v 3 ) = v 1 (α s 2 ψ + β c 2 ψ ) + v 2 c φ c ψ s ψ (α -β) + v 3 s φ c ψ s ψ (α -β) ∂ 2 H(v 1 , v 2 , v 3 ) = v 1 (α -β) c φ c ψ s ψ + v 2 (α c 2 φ c 2 ψ + β c 2 φ c 2 ψ + β s 2 φ ) +v 3 (α -β) c 2 ψ c φ s φ ∂ 3 H(v 1 , v 2 , v 3 ) = v 1 (α -β) s φ c ψ s ψ + v 2 (α -β) c 2 ψ c φ s φ + v 3 (α s 2 φ c 2 ψ +β s 2 φ s 2 ψ + β c 2 φ )
We treat three cases, at fixed u 1 , u 2 and u 3 , respectively.

Case 1. We minimize H(u 1 , •, •) on R 2 . This is equivalent to solving

∂ 2 H(u 1 , v 2 , v 3 ) = 0 ∂ 3 H(u 1 , v 2 , v 3 ) = 0.
The solution is given by

v 2 = u 1 (β -α) c φ c ψ s ψ α c 2 ψ + β s 2 ψ v 3 = u 1 (β -α) s φ c ψ s ψ α c 2 ψ + β s 2 ψ .
The minimum then reads inf w2,w3

H(u 1 , w 2 , w 3 ) = H(u 1 , v 2 , v 3 ) = αβ u 2 1 (α c 2 ψ + β s 2 ψ ) 2 (β s 2 ψ + α c 2 ψ ) = u 2 1 αβ α c 2 ψ + β s 2 ψ . (4.4)
Case 2. We minimize H(•, u 2 , •) on R 2 . This is equivalent to solving

∂ 1 H(v 1 , u 2 , v 3 ) = 0 ∂ 3 H(v 1 , u 2 , v 3 ) = 0.
The solution is given by

v 1 = u 2 (β -α) c ψ c φ s ψ c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ v 3 = u 2 (β -α) s φ c φ c 2 ψ c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ .
The minimum then reads inf w1,w3

H(w 1 , u 2 , w 3 ) = H(v 1 , u 2 , v 3 ) = αβ u 2 2 ( c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ ) 2 (β c 2 φ c 2 ψ + α s 2 ψ c 2 φ + α s 2 φ ) = u 2 2 αβ c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ . ( 4
H(w 1 , w 2 , u 3 ) = u 2 3 αβ s 2 φ (α s 2 ψ + β s 2 ψ ) + α c 2 φ . (4.6)
Slicing argument and convexity inequalities. In the three-dimensional case, the slicing argument and inequality (4.5) imply, with obvious notation,

λ 2 = (0,1) 2 α E 2 1 (y1,y3) |∇u| 2 dy 2 + β (0,1) 2 \E 2 2 (y1,y3) |∇u| 2 dy 2 + E 2 2 (y1,y3)\E 2 1 (y1,y3) \E 2 3 (y1,y3)   β 0 0 0 α 0 0 0 α   R(y)∇u, R(y)∇u dy 2 + E 2 3 (y1,y3)   α 0 0 0 β 0 0 0 β   R(y)∇u, R(y)∇u dy 2   dy 1 dy 3 ≥ (0,1) 2 inf α E 2 1 (y1,y3) v ′ (y 2 ) 2 dy 2 + β (0,1) 2 \E 2 2 (y1,y3) v ′ (y 2 ) 2 dy 2 + E 2 2 (y1,y3)\E 2 1 (y1,y3) \E 2 3 (y1,y3) v ′ (y 2 ) 2 αβ c 2 φ (β s 2 ψ + α c 2 ψ ) + β s 2 φ dy 2 + E 2 3 (y1,y3) v ′ (y 2 ) 2 αβ c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ dy 2 , v ∈ H 1 0 ((0, 1)) dy 1 dy 3
Similar bounds can be derived for λ 1 and λ 3 .

Next, we can use Proposition 2.5. In addition, Proposition 2.6 also holds in any dimension. Therefore

1 θ 1 + 1 θ 2 + 1 θ 3 ≤ 3θ 1 1 α + 3(1 -θ 2 ) 1 β + E3 α c 2 ψ + β s 2 ψ αβ + c 2 φ (α s 2 ψ + β c 2 ψ ) + α s 2 φ αβ + s 2 φ (α s 2 ψ + β s 2 ψ ) + α c 2 φ αβ dy + (E2\E1)\E3 β c 2 ψ + α s 2 ψ αβ + c 2 φ (β s 2 ψ + α c 2 ψ ) + β s 2 φ αβ + s 2 φ (β s 2 ψ + α s 2 ψ ) + β c 2 φ αβ dy = 3θ 1 1 α + 3(1 -θ 2 ) 1 β + θ 3 1 α + 2 β + (θ 2 -θ 1 -θ 3 ) 1 β + 2 α = 3 θ α + 1 -θ β ,
which is the desired harmonic bound.

Optimality of the bounds.

Using the natural extensions of the twodimensional test functions of Section 3 to three dimensions, one can show the harmonic and arithmetic bounds to be attained. The entire convex set is also attained by layering.

5.

Interpretation in terms of G-closure. In this last section, we briefly recall a more classic way to derive the arithmetic and harmonic bounds, also known as the Voigt and Reuss bounds. We then interpret our results in terms of quasiconvex envelope and attainment of infima following the seminal papers of Kohn and Strang [START_REF] Kohn | Optimal design and relaxation of variational problems, I, II and III[END_REF], based on the vectorization of the problem by Tartar [START_REF] Murat | H-convergence[END_REF][START_REF] Murat | Séminaire d'Analyse Fonctionnelle et Numérique[END_REF], and further studied by Allaire and Francfort [START_REF] Allaire | Existence of minimizers for non-quasiconvex functionals arising in optimal design[END_REF].

Derivation of the trivial bounds.

Let n ∈ N, ξ ∈ M n (R) and consider

H 1 0 ((0, 1) n , R n ) ∋ u → F (ξ, ∇u) = (0,1) n A(y)(ξ + ∇u), ξ + ∇u dy
The arithmetic bound is a consequence of inf{F (ξ, ∇u) : u ∈ H 1 0 ((0, 1) n )} ≤ F (ξ, 0) for ξ = 1 0 0 1 . The derivation of the harmonic bound is related to the same minimization on an enlarged space, namely on

L 2 # ((0, 1) n ) n×n = v ∈ L 2 ((0, 1) n ) n×n , (0,1) n v dy = 0 .
In this case, a direct computation shows inf{F (ξ, ∇u)

: u ∈ H 1 0 ((0, 1) n } ≥ inf{F (ξ, v) : v ∈ L 2 # ((0, 1) n ) n×n } = (0,1) n A -1 -1 ξ, ξ .
This inequality has an interpretation in terms of quadratic forms since the left-hand side is a quadratic form (quadratic forms are closed by Γ-convergence). Taking the inverse on both sides and reversing the inequality yields the harmonic bound for

ξ = 1 0 0 1 .
This way of deriving the trivial bounds is more straightforward and holds for any dimension. However it does not immediately suggest the form of optimal configurations, as opposed to the present 'pointwise' approach.

Interpretation in terms of quasiconvexification.

Following the approach initiated in [START_REF] Kohn | Optimal design and relaxation of variational problems, I, II and III[END_REF], we may look for the harmonic bound as the minimum of a functional depending on the geometry of the different phases of the material in (0, 1) n .

Let χ E1 and χ E2 be the characteristic functions associated to the two isotropic phases A 1 = α 0 0 α and A 2 = β 0 0 β . And let A 3 be the anisotropic phase.

For ξ ∈ M n (R), the objective is to minimize the functional

J(χ 1 , χ 2 ) = (0,1) n (χ 1 A 1 + χ 2 A 2 + (1 -χ 1 -χ 2 )A 3 ) (ξ + ∇u χ1,χ2 ), ξ + ∇u χ1,χ2 dy
on the space of admissible characteristic functions

Ξ = (χ 1 , χ 2 ) ∈ L ∞ ((0, 1) n ; {0, 1} 2 ) : χ 1 χ 2 = 0, (0,1) n χ 1 + 1 2 (1 -χ 1 -χ 2 ) dy = θ, (0,1) n χ 2 + 1 2 (1 -χ 1 -χ 2 ) dy = 1 -θ ,
θ ∈ (0, 1) being fixed. Here u χ1,χ2 realizes the minimum of

F (ξ, ∇(•)) on H 1 0 ((0, 1) n ), with A(y) = χ 1 A 1 + χ 2 A 2 + (1 -χ 1 -χ 2 )A 3 . We then have inf Ξ J(χ 1 , χ 2 ) = inf Ξ inf u∈H ξ (0,1) n (χ 1 A 1 + χ 2 A 2 + (1 -χ 1 -χ 2 )A 3 ) ∇u, ∇u dy = inf u∈H ξ (0,1) n inf χ1,χ2=0,1;χ1χ2=0 (χ 1 A 1 + χ 2 A 2 +(1 -χ 1 -χ 2 )A 3 )∇u, ∇u + l 1 (χ 1 + 1 2 (1 -χ 1 -χ 2 )) +l 2 (χ 2 + 1 2 (1 -χ 1 -χ 2 )) dy,
where

H ξ = {u ∈ H 1 ((0, 1) n , R n ) : u(x) = ξ • x on ∂(0, 1
) n } and l 1 , l 2 are Lagrange multipliers to impose the volume averages of χ 1 and χ 2 .

Therefore, one can interpret the computation of the optimal bounds on the trace of the homogenized matrix as the computation of an infimum of a vector-valued integral functional. In order to apply the direct method of the calculus of variations, one may first compute the quasiconvex envelope of the integrand. The computation of a quasiconvex envelope is rarely explicit and this is the main reason why it may seems surprising at a first glance that we achieve this difficult task with such very simple arguments.

Actually, the bounds we have obtained being the trivial bounds, the quasiconvex envelope coincides with the convex envelope and inf Ξ J(χ 1 , χ 2 ) = inf v∈V ξ (0,1) n inf χ1,χ2=0,1;χ1χ2=0

(χ 1 A 1 + χ 2 A 2 + (1 -χ 1 -χ 2 )A 3 )v, v +l 1 (χ 1 + 1 2 (1 -χ 1 -χ 2 )) + l 2 (χ 2 + 1 2 (1 -χ 1 -χ 2 )) dy,
where V ξ = {v ∈ L 2 # ((0, 1) n×n : (0,1) n v = ξ}. There exists a family of matrix fields for which this infimum is attained and which saturate the harmonic bound. The associated minimizers v ∈ V ξ turn out to belong to H ξ . The matrices are represented on Figure 3.2 with varying θ 1 and θ 2 such that θ 1 + θ 2 = 2θ.

Turning now to the arithmetic bounds, there also exists a family of matrix fields that saturate this inequality. These matrices are represented on Figure 3.1, with varying θ 1 and θ 2 such that θ 1 + θ 2 = 2θ. In particular, they are divergence free. Theorem 3. Let {α 1 , . . . , α N } be N ordered conductivities. We denote by A = {A β1,...,βn } the set of diagonal matrices of order n such that A β1,...,βn (i, i) = α βi . Let A : (0, 1) n → R n be such that, up to a rotation and for almost all x ∈ (0, 1) n , A(x) ∈ A, we set A = β1,...,βn=1,N χ β1,...,βn R(x) T A β1,...,βn R(x), where χ β1,...,βn are characteristic functions and R(x) are rotations. We denote by θ β1,...,βn = (0,1) n χ β1,...,βn . Let now fix to θ αi ∈ (0, 1) the proportions of the conductivities α i .

We have i=1,N θ αi = 1. Consider the following relationships between the phases, An example of optimal configuration is given for the arithmetic bounds in Figure 5.1. The rotation of π/2 of this configuration gives an optimal configuration for the harmonic bound.

Fig. 1 . 1 .

 11 Fig. 1.1. Comparison of the bounds for discrete and isotropic continuous materials

Fig. 1 .

 1 Fig. 1.2. a conducting network and the corresponding continuous coefficients

  Fig. 1.3. Polycrystal

Fig. 3 . 1 .

 31 Fig. 3.1. Optimal geometry for the arithmetic bound

. 5 ) 3 .

 53 Case Case 3 is the same as Case 2 up to switching c φ and s φ . Thus inf w1,w2

5. 3 .

 3 Extensions to any number of conductivities and dimensions. The interpretation in terms of quasiconvexification and the generality of the Voigt and Reuss bounds allow us to state a general version of Theorem 1.

  ...,1,βn + • • • + β2,...,βn =1 n -1 n θ 1,β2,...,βn = θ α1 . . . θ N,...,N + βn =N 1 n θ N,...,N,βn + • • • + β2,...,βn =N n -1 n θ N,β2,...,βn = θ αN , the associated G-closure is the bounded set of λ 1 , . . . , λ n satisfying simultaneously all the following inequalities in eigenvalue space 

Fig. 5 . 1 .

 51 Fig. 5.1. Optimal geometry for the arithmetic bound in 2 dimensions for 4 conductivities

3 Fig. 5 . 2 .

 352 Fig. 5.2. Exact bounds in 3 dimensions (note that 1 ≤ λ i ≤ 3)

Acknowledgements. We thank Enzo Nesi for valuable comments and help concerning the link of the present work with the theory of G-closure, and acknowledge support from the Marie Curie Research Training Network MRTN-CT-2004-505226 'Multi-scale modelling and characterisation for phase transformations in advanced materials' (MULTIMAT).


Roberto Alicandro

Marco Cicalese

Antoine Gloria

Philippe Ciarlet

Mathematical derivation of a rubber-like stored energy functional

In this note, we consider a stochastic network of interacting points to which we associate an energy. We study the variational convergence of such an energy when the typical distance of the network goes to zero. We prove that the limit energy can be written as an integral functional, whose energy density is deterministic, hyperelastic and frame-invariant. This derivation allows us in particular to obtain a continuous energy density associated to cross-linked polymer networks.

Introduction

As a starting point, we consider a particular case of the stochastic networks introduced by Blanc, Le Bris and Lions in [START_REF] Blanc | The energy of some microscopic stochastic lattices[END_REF][START_REF] Blanc | Du discret au continu pour des réseaux aléatoires d'atomes[END_REF]. In what follows, (Ω, F, P ) denotes a given probability space. Let Λ = {x i } i∈Z d ∈ R d Z d be a set of points. We say that Λ is an admissible set of points if it satisfies the three following conditions:

Email addresses: alicandr@unicas.it (Roberto Alicandro), cicalese@unina.it (Marco Cicalese), gloria@cermics.enpc.fr (Antoine Gloria).

i. there exists R > 0 such that #Λ ∩ B(x, R) > 0 for all x ∈ R d ; ii. there exist R ′ > 0 and k ∈ N such that #Λ ∩ B(x, R ′ ) ≤ k for all x ∈ R d ; iii. there exists r > 0 such that d(x i , Λ \ {x i }) ≥ r for all i ∈ Z d .

In particular, to each admissible set of points Λ one can associate a Delaunay triangulation D(Λ).

A stochastic lattice L : Ω → R d Z d is said to be admissible, if for P -almost every ω ∈ Ω, L(ω) is an admissible set of points, and if the Delaunay triangulation is regular in the sense of the interpolation theory. Given a measure preserving group {τ k } k∈Z d acting on Ω, L is said to be stationary if L(τ k (ω)) = L(ω)k for all k ∈ Z d and ergodic if, for any A ∈ F, τ k (A) = A implies P (A) = 0 or P (A) = 1. Note that in two dimensions, any admissible set of points has a regular Delaunay triangulation. In three dimensions, the conditions i., ii. and iii. are not sufficient to ensure the existence of a regular triangulation.

To each realization ω of the stochastic lattice, we associate an energy which only depends on the realization. For all regular bounded open subset A of R d , and all u : L(ω) → R n , we define the energy of the lattice by

E(u, A)(ω) = xi∈L(ω)∩A x j = x i ∈ L(ω) ∩ A [x i , x j ] ⊂ A J(x j -x i )f u(x j ) -u(x i ) |x j -x i | , (1)
where J : R d → R and f : R n → R are given functions that do not depend on ω and A and where

[x i , x j ] = tx i + (1 -t)x j , t ∈ [0, 1]
. The precise forms of J and f will be prescribed in the following section.

In what follows we consider a scaled version of the energy in (1); i.e. for all ε > 0, and

u ε : εL(ω) → R n , we set E ε (u ε , A)(ω) = xi∈εL(ω)∩A ε d x j = x i ∈ εL(ω) ∩ A [x i , x j ] ⊂ A J x j -x i ε f u ε (x j ) -u ε (x i ) |x j -x i | . (2)
2. Derivation of a macroscopic model

Main result

We make the following set of assumptions (H) on the energy E ε :

(H)

             There exist p > 1, C 1 , C 2 > 0, such that C 1 |v| p -1 ≤ f (v) ≤ C 2 |v| p + 1, ∀v ∈ R n J ≥ 0 and inf B(0,4R) J(z) > 0 R d J(z)dz < ∞
In particular, if n = d and for all w ∈ S 1 (R d ), f (w) = 0, then u : x → x is a minimizer of the discrete energy and the lattice is a ground state of the system of interacting points. This may model in particular the ground state of a rubber-like polymer at the scale of the polymer chain (see [START_REF] Treolar | The physics of rubber elasticity[END_REF]).

For all ω ∈ Ω, let T = {T i } i∈Z d be the set of the triangles of a regular Delaunay triangulation associated to εL(ω). In what follows, in order to study the asymptotics of the energy in (2) as ε → 0 in the framework of Γ-convergence (see [START_REF]Braides Γ-convergence for beginners[END_REF] for an introduction to the subject), we identify each u : εL(ω) → R n with its continuous piecewise affine interpolation on T .

Theorem 1 Let D be an open bounded subset of R d , L : Ω → R d Z d
be an admissible stationary and ergodic stochastic lattice, and let f and J satisfy Hypotheses (H). For P -almost all ω ∈ Ω, E ε (•, D)(ω) Γ -L p (D) converges to E(•, D) which is finite on W 1,p (D) and defined by

E(v, D) = D W hom (∇v(x))dx, (3)
where W hom is an homogeneous in space quasiconvex energy density satisfying a growth condition of order p, and given by the following asymptotic formula

W hom (ξ) = lim N →∞ 1 N d Ω inf E 1 (u, (0, N ) d )(ω), u(z) = ξ • z for z ∈ L(ω), dist(z, ∂(0, N ) d ) ≤ 2R dP (ω), (4)
In addition, infimum problems with prescribed boundary conditions also converge.

Sketch of the proof

The details of the following proof will be given in [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF].

Individual compactness. The first step of the proof consists in obtaining individual compactness. For all ω ∈ Ω, we prove that there exists a subsequence ε n , such that

E εn (•, D)(ω) Γ-converges to some E(•, D)(ω) on W 1,p (D)
, for any open bounded regular subset D ⊂ R d . To this end, one proceeds as in [START_REF] Alicandro | A general integral representation result for the continuum limits of discrete energies with superlinear growth[END_REF].

The estimates are however more delicate to get. Fixed a regular Delaunay triangulation associated to ω, one first singles out the contribution to the energy due to the interactions between nearest-neighbors, by which we mean two vertices of the same triangle, and, using the coercivity hypothesis on f , obtain the coercivity of the Γ-liminf. Up to introducing a coarser periodic reference lattice, one can classify the interactions according to their range. This allows us to adapt the argument used in [START_REF] Alicandro | A general integral representation result for the continuum limits of discrete energies with superlinear growth[END_REF] where the periodic-lattice case is considered in the deterministic setting and to get the growth condition from above and the measure property as set function of the Γ-limsup E ε (u, A) for any regular open set A ⊂ R d . Using then the compactness property of Γ-convergence and the integral representation result in [START_REF] Buttazzo | Integral representation and relaxation of local functionals[END_REF], we prove that, up to extractions, there exists a Γ-limit E(•, D)(ω), which can be expressed in terms of an integral functional. The second part of the proof is devoted to show that E(•, D)(ω) does not depend on ω and that the associated energy density W hom does not depend on the space variable.

Ergodicity of subadditive processes. The result is achieved if we prove the existence of the limit (4) for the asymptotic formula. The proof relies on the subadditive ergodic theorem, and more precisely on the variant used in [START_REF] Maso | Nonlinear stochastic homogenization and ergodic theory[END_REF]. The energy E 1 (•, D) is not subadditive since the interactions are not local. So is its infimum. One may however introduce a modified energy defined on the borelians of R d , which is subadditive and such that the asymptotic limits of type (4) coincide for both energies, the convergence of one implying the convergence of the other. Doing so, the ergodic theorem shows the existence of ( 4) and therefore concludes the proof. Let us quickly define the modified energy in the case of finite-range interactions, namely ∃L ∈ N :

J(z) = 0 if |z| > L. Let B ∈ B(R d
G ξ (B)(ω) = inf{E(u, C(B))(ω), u(x i ) = ξ • x if d(x i , ∂C(B)) ≤ R} + α perim(C(B)),
where perim denotes the perimeter of the set. For all ξ, there exists α big enough such that G ξ is a subadditive set functional. In addition, for regular sets B, lim G ξ (tB)(ω) and its independence upon ω, which concludes the proof.

Remarks and extensions

Some comments are in order. First, it is easily shown that W hom is frame-invariant. Second, if there exists a stationary action group of rotations, then W hom is isotropic.

Volume effects are not taken into account in Theorem 1. Given a regular Delaunay triangulation associated to the lattice, one can slightly modify the energy of the lattice by adding a term penalizing volume changes on each triangle of T and satisfying a growth condition of order p from above. In such a case, Theorem 1 also holds. However, this modification does not cover incompressible behaviours and does not preserve the orientation. To obtain more realistic energy densities, one may proceed in two steps: first use a cut-off η on the 'realistic' volumetric energy (see [7]) in order to make it satisfy a suitable growth condition and pass to the limit as ε → 0, and then let the cut-off parameter η go to infinity to recover the 'original' volumetric properties. However it is not clear whether it is possible to switch the limits. These remarks and extensions will be detailed in a forthcoming paper [START_REF] Alicandro | Integral representation results for energies defined on stochastic lattices and application to nonlinear elasticity[END_REF].

We finally point out that the variational convergence of stochastic energies defined on a one-dimensional periodic network of points in interaction has been considered in [START_REF] Iosifescu | Variational limit of a one-dimensional discrete and statistically homogeneous system of material points[END_REF][START_REF] Iosifescu | Variational limit of a one-dimensional discrete and statistically homogeneous system of material points[END_REF].

  ) be a bounded borelian set. We denote by C(B) the union of the finite number of closed unit cubes C centered in points of Z d such that C ∩ B = ∅, and by C(B) the interior of C(B). For all α > 0 and for all ξ ∈ M d (R), we set


Evan Senter
email: evansenter@gmail.com

Saad Sheikh
email: sheikh@cise.ufl.edu

Ivan Dotu
email: ivan.dotu@bc.edu

Yann Ponty
email: yann.ponty@lix.polytechnique.fr

Peter Clote
email: clote@bc.edu

Using the Fast Fourier Transform to accelerate the computational search for RNA conformational switches (extended abstract)

Keywords: RNA secondary structure, partition function, fast Fourier transform, Lagrange interpolation

We describe the broad outline of a new thermodynamicsbased algorithm, FFTbor, that uses the fast Fourier transform to perform polynomial interpolation to compute the Boltzmann probability that secondary structures differ by k base pairs from an arbitrary reference structure of a given RNA sequence. The algorithm, which runs in quartic time O(n 4 ) and quadratic space O(n 2 ), is used to determine the correlation between kinetic folding speed and the ruggedness of the energy landscape, and to predict the location of riboswitch expression platform candidates.

In [START_REF] Freyhult | Boltzmann probability of RNA structural neighbors and riboswitch detection[END_REF], we developed a dynamic programming algorithm, RNAbor, which simultaneously computes for each integer k, the Boltzmann probability p k = Z k Z of the subensemble of structures whose base pair distance to a given initial, or reference, structure S * is k. 4 RNAbor stores the value of the (partial) partition functions Z k (i, j) for all 1 ≤ i ≤ j ≤ n and 0 ≤ k ≤ n, each of which requires quadratic time to compute. Thus it follows that RNAbor runs in time O(n 5 ) and space O(n 3 ), which severely limits its applicability to genomic annotation. This restriction is somewhat mitigated by the fact that in [START_REF] Clote | Maximum expected accuracy structural neighbors of an RNA secondary structure[END_REF], we showed how to use sampling to efficiently approximate RNAbor in cubic time O(n 3 ) and quadratic space O(n 2 ), provided that the starting structure S * is the minimum free energy (MFE) structure. We expect that a more efficient version of RNAbor could be used in applications in genomics and synthetic biology, to detect potential conformational switches -RNA sequences containing two or more (distinct) metastable structures. is taken to be the minimum free energy structure. (Bottom) Expected base pair distance k k • Z k /Z from the reference structure of the guanine riboswitch of Bacillus subtilis, depicted in Figure 1A of [START_REF] Serganov | Structural basis for discriminative regulation of gene expression by adenine-and guanine-sensing mRNAs[END_REF]. FFTbor was run on all 144 nt windows of CP003901.1/847300-848050, comprising the 5 ′ untranslated region of the XPT gene (guanosine monophosphate reductase, with coding region at CP003901.1/848026-848607) of the unrelated organism Streptococcus pyogenes A20. FFTbor detects the guanine riboswitch at position 847848, where expected base pair distance to S * is minimized (53.79) corresponding to a Z-score of -6.95. This prediction corresponds well with the Rfam prediction at nearby position 847844.

In this abstract, we announce a radically different algorithm, FFTbor, that uses polynomial interpolation to compute the coefficients p 0 , . . . , p n-1 of the polynomial p(x) = p 0 +p 1 x+• • •+p n-1 x n-1 , where p k is defined by p k = Z k Z . Due to severe numerical instability issues in both the Lagrange interpolation formula and in Gaussian elimination, we employ the Fast Fourier Transform (FFT) to compute the inverse Discrete Fourier Transform (DFT) on values y 0 , . . . , y n-1 , where y k = p(ω k ) and ω = e 2πi/n is the principal nth complex root of unity. This gives rise to an improved version of RNAbor, denoted FFTbor, which runs in time O(n 4 ) and space O(n 2 ) on a single processor, and in time O(n 3 ) on a theoretical n-cores processor or cluster (e.g. using OpenMP). Figure 1 (top) depicts the rugged energy landscape typical of a conformational switch, while Figure 1 (bottom) depicts expected base pair distance, of each size 144 window in the 5 ′ -UTR of S. pyogenes, to the XPT riboswitch structure of B. subtilis.

Fig. 1 .

 1 Fig. 1. (Top) Output of FFTbor on the 46 nt spliced leader conformational switch of Leptomonas collosoma, where reference structure S * is taken to be the minimum free energy structure. (Bottom) Expected base pair distancek k • Z k /Z from the reference structure of the guanine riboswitch of Bacillus subtilis, depicted in Figure1Aof[START_REF] Serganov | Structural basis for discriminative regulation of gene expression by adenine-and guanine-sensing mRNAs[END_REF]. FFTbor was run on all 144 nt windows of CP003901.1/847300-848050, comprising the 5 ′ untranslated region of the XPT gene (guanosine monophosphate reductase, with coding region at CP003901.1/848026-848607) of the unrelated organism Streptococcus pyogenes A20. FFTbor detects the guanine riboswitch at position 847848, where expected base pair distance to S * is minimized (53.79) corresponding to a Z-score of -6.95. This prediction corresponds well with the Rfam prediction at nearby position 847844.

Here Z denotes the partition function, defined as the sum of all Boltzmann factors exp(-E(S)/RT ), over all secondary structures S of a given RNA sequence, R denotes the universal gas constant and T absolute temperature. Similarly Z k denotes the sum of all Boltzmann factors of all structures S, whose base pair distance to the initial structure S * is exactly k.


Vladimir Reinharz

Yann Ponty

Jérôme Waldispühl

A linear inside-outside algorithm for correcting sequencing errors in structured RNA sequences

Keywords: RNA, mutations, secondary structure

Analysis of the sequence-structure relationship in RNA molecules are essential to evolutionary studies but also to concrete applications such as error-correction methodologies in sequencing technologies. The prohibitive sizes of the mutational and conformational landscapes combined with the volume of data to proceed require efficient algorithms to compute sequence-structure properties. More specifically, here we aim to calculate which mutations increase the most the likelihood of a sequence to a given structure and RNA family. In this paper, we introduce RNApyro, an efficient linear-time and space inside-outside algorithm that computes exact mutational probabilities under secondary structure and evolutionary constraints given as a multiple sequence alignment with a consensus structure. We develop a scoring scheme combining classical stacking base pair energies to novel isostericity scales, and apply our techniques to correct point-wise errors in 5s rRNA sequences. Our results suggest that RNApyro is a promising algorithm to complement existing tools in the NGS error-correction pipeline.

Introduction

Ribonucleic acids (RNAs) are found in every living organisms and exhibit a broad range of functions, from catalyzing chemical reactions as the RNase P or the group II introns, hybridizing messenger RNA to regulate gene expression, to ribosomal RNA (rRNA) synthesizing proteins. Those functions require specific structures, encoded in their nucleotide sequence. Although the functions, and thus the structures, need to be preserved through various organisms, the sequences can greatly differ from one organism to another. This sequence diversity coupled with the structural conservation is a fundamental asset for evolutionary studies. To this end, algorithms to analyze the relationship between RNA mutants and structures are required.

For half a century, biological molecules have been studied as a proxy to understand evolution [START_REF] Zuckerkandl | Molecules as documents of evolutionary history[END_REF], and due to their fundamental functions and remarkably conserved structures, rRNAs have always been a prime candidate for phylogenetic studies [START_REF] Olsen | Microbial ecology and evolution: a ribosomal RNA approach[END_REF][START_REF] Olsen | Ribosomal RNA: a key to phylogeny[END_REF]. In recent years, studies as the Human Microbiome Project [START_REF] Peter J Turnbaugh | The Human Microbiome Project[END_REF] benefited of new technologies such as the NGS techniques to sequence as many new organisms as possible and extract an unprecedented flow of new information. Nonetheless, these high-throughput techniques typically have high error rates that make their applications to metagenomics (a.k.a. environmental genomics) studies challenging. For instance, pyrosequencing as implemented by Roche's 454 produces may have a error rate raising up to 10%. Because there is no cloning step, resequencing to increase accuracy is not possible and it is therefore vital to disentangle noise from true sequence diversity in this type of data [START_REF] Quince | Accurate determination of microbial diversity from 454 pyrosequencing data[END_REF]. Errors can be significantly reduced when large multiple sequence alignments with close homologs are available, but in studies of new or not well known organisms, such information is rather sparse. In particular, it is common to is not enough similarity to differentiate between the sequencing errors and the natural polymorphisms that we want to observe, often inflating the diversity estimates [START_REF] Victor Kunin | Wrinkles in the rare biosphere: pyrosequencing errors can lead to artificial inflation of diversity estimates[END_REF]. A few techniques have been developed to remedy to this problem [START_REF] Ar Quinlan | Pyrobayes: an improved base caller for SNP discovery in pyrosequences[END_REF][START_REF] Medvedev | Error correction of highthroughput sequencing datasets with non-uniform coverage[END_REF] but they do not take into account all the available information. It is therefore essential to develop methods that can exploit any type of signal available to correct errors.

In this paper, we introduce RNApyro, a novel algorithm to that enable us to calculate precisely mutational probabilities in RNA sequences with a conserved consensus secondary structure. We show how our techniques can exploit the structural information embedded in physics-based energy models, covariance models and isostericity scales to identify and correct point-wise errors in RNA molecules with conserved secondary structure. In particular, we hypothesize that conserved consensus secondary structures combined with sequence profiles provide an information that allow us to identify and fix sequencing errors.

Here, we expand the range of algorithmic techniques previously introduced with RNAmutants [START_REF] Waldispühl | Efficient Algorithms for Probing the RNA Mutation Landscape[END_REF]. Instead of exploring the full conformational landscape and sample mutants, we develop an insideoutside algorithm that enables us to explore the complete mutational landscape with a fixed secondary structure and to calculate exactly mutational probability values. In addition to a gain into the numerical precision, this strategy allows us to drastically reduce the computational complexity (O(n 3 • m 2 ) for the original version of RNAmutants to O(n • m 2 ) for RNApyro, where n is the size of the sequence and m the number of mutations).

We design a new scoring scheme combining nearest-neighbor models [START_REF] Douglas | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF] to isostericity metrics [START_REF] Stombaugh | Frequency and isostericity of RNA base pairs[END_REF]. Classical approaches use a Boltzmann distribution whose weights are estimated using a nearest-neighbour energy model [START_REF] Douglas | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]. However, the latter only accounts for canonical and wobble, base pairs. As was shown by Leontis and Westhof [START_REF] Leontis | Geometric nomenclature and classification of RNA base pairs[END_REF], the diversity of base pairs observed in tertiary structures is much larger, albeit their energetic contribution remains unknown. To quantify geometrical discrepancies, an isostericity distance has been designed [START_REF] Stombaugh | Frequency and isostericity of RNA base pairs[END_REF], increasing as two base pairs geometrically differ from each other in space. Therefore, we incorporate these scores in the Boltzmann weights used by RNApyro.

We illustrate and benchmark our techniques for point-wise error corrections on the 5S ribosomal RNA. We choose the latter since it has been extensively used for phylogenetic reconstructions [START_REF] Hori | Origin and Evolution of Organisms as Deduced from 5s Ribosomal RNA Sequences[END_REF] and its sequence has been recovered for over 712 species (in the Rfam seed alignment with id RF00001). Using a leave one out strategy, we perform random distributed mutations on a sequence. While our methodology is restricted to the correction of point-wise error in structured regions (i.e. with base pairs), we show that RNApyro can successfully extract a signal that can be used to reconstruct the original sequence with an excellent accuracy. This suggests that RNApyro is a promising algorithm to complement existing tools in the NGS error-correction pipeline.

The algorithm and the scoring scheme are presented in Sec. 2. Details of the implementation and benchmarks are in Sec. 3. Finally, we discuss future developments and applications in Sec. 4.

Methods

We introduce a probabilistic model, which aims at capturing both the stability of the folded RNA and its ability to adopt a predefined 3D conformation. To that purpose, a Boltzmann weighted distribution is used, based on a pseudo-energy function E(•) which includes contributions for both the free-energy and its putative isostericity towards a multiple sequence alignment. In this model, the probability that the nucleotide at a given position needs to be mutated (i.e. corresponds to a sequencing error) can be computed using a variant of the Inside-Outside algorithm [START_REF] Lari | The estimation of stochastic context-free grammars using the Inside-Outside algorithm[END_REF].

Probabilistic model

Let Ω be an gap-free RNA alignment sequence, S its associated secondary structure, then any sequence s is assigned a probability proportional to its Boltzmann factor

B(s) = e -E(s) RT , with E(s) := α • ES(s, S) + (1 -α) • EI(s, S, Ω),
where R is the Boltzmann constant, T the temperature in Kelvin, ES(s) and EI(s, S, Ω) are the free-energy and isostericity contributions respectively (further described below), and α ∈ [0, 1] is an arbitrary parameter that sets the relative weight for both contributions.

Energy contribution

The free-energy contribution in our pseudo-energy model corresponds to an additive stacking-pairs model, taking values from the Turner 2004 model retrieved from the NNDB [START_REF] Douglas | NNDB: the nearest neighbor parameter database for predicting stability of nucleic acid secondary structure[END_REF]. Given a candidate sequence s for a secondary structure S, the free-energy of S on s is given by

ES(s, S) = (i,j)→(i ,j )∈S stacking pairs ES β s i s j →s i s j
where ES β ab→a b is set to 0 if ab = ∅ (no base-pair to stack onto), the tabulated free-energy of stacking pairs (ab)/(a b ) in the Turner model if available, or β ∈ [0, ∞] for non-Watson-Crick/Wobble entries (i.e. not in {GU, UG, CG, GC, AU or UA}). This latter parameter allows to choose whether to simply penalize invalid base pairs, or forbid them altogether (β = ∞). The loss of precision due to this simplification of the Turner model remains reasonable since the targeted secondary structure is fixed (e.g. multiloops do not account for base-specific contributions). Furthermore, it greatly eases the design and implementation of dynamic-programming equations.

Isostericity contribution

The concept of isostericity score [START_REF] Stombaugh | Frequency and isostericity of RNA base pairs[END_REF] is based on the geometric discrepancy (superimposability) of two base-pairs, using individual additive contributions computed by Stombaugh et al [START_REF] Stombaugh | Frequency and isostericity of RNA base pairs[END_REF]. Let s be a candidate sequence for a secondary structure S, given in the context of a gap-free RNA alignment Ω, we define the contribution of the isostericity to the pseudo-energy as ES(s, S, Ω) = (i,j)∈S pairs

EI Ω (i,j),s i s j , where EI Ω (i,j),ab :=

s ∈Ω ISO((s i , s j ), (a, b)) |Ω|
is the average isostericity of a base-pair in the candidate sequence, compared with the reference alignment. The ISO function uses the Watson-Crick/Watson-Crick cis isostericity matrix computed by Stombaugh et al [START_REF] Stombaugh | Frequency and isostericity of RNA base pairs[END_REF]. Isostericity scores range between 0 and 9.7, with 0 being assigned to a perfect isostericity, and a penalty of 10 is used for missing entries. The isostericity contribution will favor exponentially sequences that are likely to adopt a similar local conformation as the sequences contained in the alignment.

Mutational profile of sequences

Let s be an RNA sequence, S a reference structure, and m ≥ 0 a desired number of mutations. We are interested in the probability that a given position contains a specific nucleotide, over all sequences having at most m mutations from s (formally

P(s i = x | s, Ω, S, m)).
We define a variant of the Inside-Outside algorithm [START_REF] Lari | The estimation of stochastic context-free grammars using the Inside-Outside algorithm[END_REF], allowing us to obtain the desired probability, the two functions Z * * and Y * * . The former, defined in Equations ( 2) and (3), is analogous to the inside algorithm. It is the partition function, i.e. the sum of Boltzmann factors, over all sequences within [i, j], knowing that position i -1 is composed of nucleotide a (resp. j + 1 is b), within m mutations of s. The latter, defined by Equations ( 4) and ( 5), computes the outside algorithm, i.e. the partition function over sequences within m mutations of s, restricted to two intervals [0, i] ∪ [j, n -1], and knowing that position i + 1 is composed of a (resp. j -1 is b). A suitable combination of these terms, given in Equation [START_REF] Medvedev | Error correction of highthroughput sequencing datasets with non-uniform coverage[END_REF], gives the total weight, and in turn the probability, of seeing a specific base at a given position.

Definitions

Let B := {A, C, G, U} be the set of nucleotides. Given s ∈ B n an RNA sequence, let s i be the nucleotide at position i. Let Ω be a set of un-gapped RNA sequences of length n, and S a secondary structure without pseudoknots. Formally, if (i, j) and (k, l) are base pairs in S, there is no overlapping extremities {i, j} ∩ {k, l} = ∅ and either the intersection is empty (

[i, j] ∩ [k, l] = ∅) or one is included in the other ([k, l] ⊂ [i, j] or [i, j] ⊂ [k, l]).
Let us then remind a Hamming distance function δ : B * × B * → N + , which takes two sequences s and s as input, |s | = |s |, and returns the number of differing positions. Finally, let us denote by E Ω,β (i,j),ab→ab the local contribution of a base-pair (i, j) to the pseudo-energy, such that

E Ω,β (i,j),ab→a b = α • ES β ab→a b + (1 -α) • EI Ω (i,j),a b . (1)
s i s j a b m i j a a b m -δ a ,si i + 1 j [Position i unpaired] a b a b [Position i paired anywhere else] m -m -δ a b ,sisk m i + 1 j k k -1 k + 1 a b a b m-m -δ a b ,sisj i + 1 j -1 [Paired ends + Stacking pairs]
Figure 1: Principle of the inside computation (aka partition function). Any sequence with m mutations over an interval [i, j] can be decomposed as a sequence over [i + 1, j] preceded by a, possibly mutated, base at i (Unpaired case), a sequence over [i + 1, j -1] surrounded by some base-pair (Stacking-pair case), or as two sequences over [i + 1, k -1] and [k + 1, j], completed by some base-pair (General base-pairing case). In each case, one has to investigate any possible ways to distribute mutations between the different sequences and locally instanciated bases.

Inside

The Inside function Z m

(i,j) [a,b]
is the partition function, i.e. the sum of Boltzmann factors, over all sequences in the interval [i, j], at distance m of s [i,j] , and having flanking nucleotides a and b (at positions i -1 and j + 1 respectively). Such terms can be defined by recurrence, for which the following initial conditions holds:

∀i ∈ [0, n -1] : Z m (i+1,i) [a,b] = 1 If m = 0 0 Otherwise. (2)
In other words, the set of sequences at distance m of the empty sequence is either empty if m > 0, or restricted to the empty sequence, having energy 0, if m = 0. Since the energetic terms only depend on base pairs, they are not involved in the initial conditions. The main recursion is composed of four terms:

Z m (i,j) [a,b] :=                                      a ∈B, δ a ,s i ≤m Z m-δ a ,s i (i+1,j) [a ,b] If S i = -1 a ,b ∈B 2 , δ a b ,s i s j ≤m e -E Ω,β (i,j),ab→a b RT • Z m-δ a b ,s i s j (i+1,j-1) [a ,b ] Elif S i = j ∧ S i-1 = j + 1 a ,b ∈B 2 , δ a b ,s i s k ≤m m-δ a b ,s i s k m =0 e -E Ω,β (i,k),∅→a b RT • Z m-δ a b ,s i s k -m (i+1,k-1) [a ,b ] • Z m (k+1,j) [b ,b] Elif S i = k ∧ i < k ≤ j 0 Otherwise (3) si sj a b m 0 i j n -1 a a b m -δ a ,si 0 i -1 j n -1 [Position i unpaired] a b a b m -δ a b ,sisj 0 i -1 j + 1 n -1 [Paired ends + Stacking pairs] a b b m m -δ a b ,sisk -m [Position i paired to the right] 0 i -1 j k -1 k k + 1 n -1 b a b [Position i paired to the left] m m -δ a b ,sisk -m 0 k -1 k k + 1 i -1 i j n -1
Figure 2: Principle of the outside computation. Note that the outside algorithm uses intermediate results from the inside algorithm, therefore its efficient implementation requires an implementation of the inside computation.

The cases can be broken down as follows:

S i = -1:
If the nucleotide at position i is unpaired, then any sequence consists in a, possibly mutated, nucleotide a at position i, followed by a sequence over [i + 1, j] having either m -δ a ,s i , accounting for a possible mutation at position i, and having flanking nucleotides a and b.

S i = j and S i-1 = j + 1: Any sequence generated in [i, j] consists of two, possibly mutated, nucleotides a and b , flanking a sequence over [i + 1, j -1] having distance m -δ a b ,s i s j (to avoid exceeding the targeted distance m). Since positions i and i -1 are paired with j and j + 1 respectively, then a stacking energy contribution is added.

S i = k and i < k ≤ j: If position i is paired and not involved in a stacking, then the only term contributing directly to the energy is the isostericity of the base pair (i, k). Any sequence on [i, j] consists of two nucleotide a and b at positions i and k respectively, flanking a sequence over interval [i + 1, k -1] and preceding a (possibly empty) sequence interval [k + 1, j]. Since the number of mutations sum to m over the whole sequence must , then a parameter m is introduced to distribute the remaining mutations between the two sequences.

Else: In any other case, we are in a derivation of the SCFG that does not correspond to the secondary structure S, and we return 0.

Outside

The Outside function, Y, is the partition function considering only the contributions of subsequences [0, i] ∪ [j, n -1] over the mutants of s having exactly m mutations between [0, i] ∪ [j, n -1] and whose nucleotide at position i + 1 is a (resp. in position j -1 it is b). We define function Y m

(i,j) [a,b]
as a recurrence, and will use as initial conditions:

Y m (-1,j) [X,X] := Z m (j,n-1) [X,X] (4)
The recurrence, as shown below, will increase the interval [i, j] by decreasing i when it is not base paired. If it is with a position k > j, we increase j to include it. Thus, when we need to evaluate an interval as (-1, j), all stems between (0, j) are taken into account and the structure between (j, n -1) must be a set of independent stems. Therefore, all the outside energy between [j, n

-1] is equal to Z m (j,n-1) [X,X]
, for any X ∈ B. The recursion itself is as follows.

Y m (i,j) [a,b] =                                                      a ∈B, δ a ,s i ≤m Y m-δ a ,s i (i-1,j) [a ,b] Elif S i = -1 a b ∈B 2 , δ a b ,s i s j ≤m e -E Ω,β (i,j),ab→a b RT • Y m-δ a b ,s i s j (i-1,j+1) [a ,b ] Elif S i = j ∧ S i+1 = j -1 a b ∈B 2 , δ a b ,s i s k ≤m m-δ a b ,s i s k m =0 e -E Ω,β (i,k),∅→a b RT • Y m-δ a b ,s i s k -m (i-1,k+1) [a ,b ] • Z m (j,k-1) [b,b ] Elif S i = k ≥ j a b ∈B 2 , δ a b ,s k s i ≤m m-δ a b ,s k s i m =0 e -E Ω,β (k,i),∅→a b RT • Y m-δ a b ,s k s i -m (k-1,j) [a ,b] • Z m (k+1,i-1) [a ,b ] Elif -1 < S i = k < i 0 Otherwise (5)
The five cases can be broked down as follows.

S i = -1: If the nucleotide at position i is not paired, then the value is the same as if we decrease the lower interval bound by 1 (i.e. i -1), and consider all possible nucleotides a at position i, correcting the number of mutants in function of δ a ,s i .

S i = j and S i+1 = j -1: If nucleotide i is paired with j and nucleotide i + 1 is paired with j -11, we are in the only case were stacked base pairs can occur. We thus add the energy of the stacking and of the isostericity of the base pair (i, j). What is left to compute is the outside value for the interval [i -1, j + 1] over all possible nucleotides a , b ∈ B 2 at positions i and j respectively.

S i = k ≥ j: If nucleotide i is paired with position k ≥ j, and is not stacked inside, the only term contributing directly to the energy is the isostericity of the base pair (i, k). We can then consider the outside interval [i -1, k + 1] by multiplying it by the the inside value of the newly included interval (i.e. [j, k -1]), for all possible values a , b ∈ B 2 for nucleotides at positions i and k respectively.

-1 < S i < i: As above but if position i is to pairing with a lower value.

Else: In all other cases, we are in a derivation of the SCFG that does not correspond to the secondary structure S, and we return 0.

Combining Inside and Outside values into point-wise mutations probabilities

By construction, the partition function over all sequences at exactly m mutations of s can be described in function of the inside term as Z m (0,n-1) [X,X]

, for any nucleotide X ∈ B or in function of the outside term, for any unpaired position k:

Z m (0,n-1) [X,X] ≡ a∈B, δ a,s[k] ≤m Y m-δ a,s[k] (k-1,k+1) [a,a]
We are now left to compute the probability that a given position is a given nucleotide. We leverage the Inside-Outside construction to immediately obtain the following 3 cases. Given i ∈ [0, n -1], x ∈ B, and M ≥ 0 a bound on the number of allowed mutations.

P(s i = x | s, Ω, S, M ) := W(i, x, s, Ω, S, M ) M m=0 Z m (0,n-1) [X,X] (6)
W(i, x, s, Ω, S, M ) =                              M m=0 Y m-δx,s i (i-1,i+1) [x,x] If S i = -1 M m=0 b∈B δ xb,s i s k ≤m m-δ xb,s i s k m =0 e -E Ω,β (i,k),∅→xb RT • Y m-δ xb,s i s k -m (i-1,k+1) [x,b] • Z m (i+1,k-1) [x,b]
If S i = k > i M m=0 b∈B δ bx,s k s i ≤m m-δ bx,s k s i m =0 e -E Ω,β (k,i),∅→bx RT • Y m-δ bx,s k s i -m (k-1,i+1) [b,x] • Z m (k+1,i-1) [b,x] If S i = k < i (7)
In every case, the denominator is the sum of the partitions function of exactly m mutations, for m smaller or equal to our target M . The numerators are divided in the following three cases.

S i = -1: If the nucleotide at position i is not paired, we are concerned by the weights over all sequences which have at position i nucleotide x, which is exactly the sum of the values of

Y m-δx,s i (i-1,i+1) [x,x]
, for all m between 0 and M .

S i = k > i:
Since we need to respect the derivation of the secondary structure S, if position i is paired, we must consider the two partition functions. The outside of the base pair, and the inside, for all possible values for the nucleotide at position k, and all possible distribution of the mutant positions between the inside and outside of the base pair. We also add the term of isostericity for this specific base pair.

S i = k < i: Same as above, but with position i pairing with a lower position.

Complexity considerations

Equations ( 3) and ( 5) can be computed using dynamic programming. (not described) of the isostericity contributions in Θ(n • |Ω|). Computing any given term therefore requires Θ(m) operations.

In principle, Θ(m • n 2 ) terms, identified by (m, i, j) triplets, should be computed. However, a close inspection of the recurrences reveals that the computation can be safely restricted to a subset of intervals (i, j). For instance, the inside algorithm only requires computing intervals [i, j] that do not break any base-pair, and whose next position j + 1 is either past the end of the sequence, or is base-paired prior to i. Similar constraints hold for the outside computation, resulting in a drastic limitation of the combinatorics of required computations, dropping from Θ(n 2 ) to Θ(n) the number of terms that need to be computed and stored. Consequently the overall complexity of the algorithm is Θ(n • (|Ω| + m 2 )) arithmetic operations and Θ(n • (|Ω| + m)) memory.

Results

Implementation

The software was implemented in Python2.7 using the mpmath [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.14)[END_REF] library for arbitrary floating point precision. The source code is freely available at https://github.com/McGill-CSB/RNApyro.

The time benchmarks were done on a MacMini 2010, 2.3GHz dual-core Intel Core i5, 8GB of RAM. Since applications of RNApyro implies a need for efficiency and scalability, we present in Table 1 typical runtimes required to compute the probabilities for every nucleotide at every positions for a vast set of parameters. For those tests, both the sequences and the target secondary structure were generated at random.

Error correction in 5s rRNA

To illustrate the potential of our algorithm, we applied our techniques to identify and correct point-wise errors in RNA sequences with conserved secondary structures. More precisely, we used RNApyro to reconstruct 5s rRNA sequences with randomly distributed mutations. This experiment has been designed to suggest further applications to error-corrections in pyrosequencing data.

We built our data set from the 5S rRNA multiple sequence alignment (MSA) available in the Rfam Database 11.0 (Rfam id: RF00001). Since our software does not currently implement gaps (mainly because scoring indels is a challenging issue that cannot be fully addressed in this work), we clustered together the sequences with identical gap locations. From the 54 MSAs without gap produced, we selected the biggest MSA which contains 130 sequences (out of 712 in the original Rfam MSA). Then, in order to avoid overfitting, we used cd-hit [START_REF] Li | Cd-hit: a fast program for clustering and comparing large sets of protein or nucleotide sequences[END_REF] to remove sequences with more than 80% of sequence similarity. This operation resulted in a data set of 45 sequences. We designed our benchmark using a leave-one-out strategy. We randomly picked one sequence from our data set and performed 12 random mutations. Our sequences have 119 nucleotides, thus the number of mutations corresponds to an error-rate of 10%. We repeated this operation 10 times. The value of β has been set to 1.5 (larger values gave similar results). To estimate the impact of the distribution of the weights between the energy term and isostericity score, we used 4 different values of α = 0, 0.5, 0.8, 1.0. Similarly, we also investigated the impact of an under-and over-estimate of the number of errors, and use values equal to 50% (6 mutations) and 200% (24 mutations) of the exact number of errors (i.e. 12).

To evaluate our method, we computed a ROC curve representing the performance of a classifier based on the mutational probabilities computed by RNApyro. More specifically, we fixed a threshold λ ∈ [0, 1], and predicted an error at position i in sequence ω if and only if the probability P (i, n) of a nucleotide n ∈ {A, C, G, U } exceeds this threshold. To correct the errors we used the set of nucleotides with a probability higher than this threshold is C(i) = {n | n ∈ {A, C, G, U } and P (i, n) > λ and n = ω[i]}, where ω[i] is the nucleotide at position i in the input sequence. We note that for the lowest thresholds multiple nucleotides can be available in C(i) to correct the sequence. Here, we remind that our aim is to estimate the potential of error-correction of RNApyro and not to develop an error-correction pipe-line, which will be the subject of further studies. Finally, we progressively varied λ between 0 and 1 to calculate the ROC curve and the area under the curve (AUC). We report our results in Figure 3.

Our data demonstrate that our algorithm exhibits interesting performance for error-correction applications. First, the AUC values (up to 0.86) indicate that a signal has been successfully extracted. This result has been achieved with errors in loop regions (i.e. without base pairing information) and thus suggests that correction rates in structured regions (i.e. base paired regions) could be even higher. Next, the optimal values of α tend to be close to 0.0. This finding suggests that at this point the information issued from the stacking energy is currently modest. However, specific examples showed improved performance using this energy term. Further studies must be conducted to understand how to make the best use it. Finally, our algorithm seems robust to the number of mutations considered. Indeed, good AUC values are achieved with low estimates of the number of errors in the sequence (c.f. 50% of the errors in Fig. 3a) and with large values (c.f. 200% of the errors in Fig. 3c). It is worth noting that scoring scheme with larger weight on the isostericity

α(α(αFigure 3 :

 3 Figure3: Performance of error-correction. Subfigures accuracy with under-estimated error rate (6 mutations), exact estimates (12 mutations) and over estimates (24 mutations). We also analyze the impact of the parameter α distributing the weights of stacking pair energies vs isostericity scores and use values ranging of α = {0, 0.5, 0.8, 1.0}. The AUC is indicated in the legend of the figures. Each individual ROC curve represent the average performance over the 10 experiments.

Table 1 :

 1 Time to compute all probabilities. The first column indicates the length and the column indexes indicate the number of mutations. α is set at 0.5, β to 1.5 and |Ω| = 44.

	|s|	Number of mutations
		6	12	24
	100	35s	238s 1023s
	300 135s 594s 2460s
		25		50
	500 5400s		21003s

Namely, the Z * * and Y * * terms are computed starting from smaller values of m and interval lengths, memorizing the results as they become available to ensure constant-time access during later stages of the computation. Furthermore, energy terms E(•) can be accessed in constant time thanks to a simple precomputation

scores (low α values) seem more robust to under-and over-estimate of the number of errors.

Conclusion

In this article we presented a new and efficient way of exploring the mutational landscape of an RNA under structural constraints, and apply our techniques to identify and fix sequencing errors. In addition, we introduce a new scoring scheme combining the nearest-neighbour energy model to new isostericity matrices in order to account for geometrical discrepancies occurring during base pair replacements. The algorithm runs in Θ(n • (|Ω| + m 2 )) time and Θ(n • (|Ω| + m)) memory, where n is the length of the RNA, m the number of mutations and Ω the size of the multiple sequence alignment.

By combining into RNApyro these two approaches, the mutational landscape exploration and the pseudo energy model, we created a tool predicting the positions yielding point-wise sequencing error and correcting them. We validated our model with the 5s rRNA, as presented in Sec. [START_REF] Victor Kunin | Wrinkles in the rare biosphere: pyrosequencing errors can lead to artificial inflation of diversity estimates[END_REF]. We observed that the models with larger weights on the isostericity seems to hold a higher accuracy on the estimation of errors. This indicates that an extractable signal is contained in the isostericity. Importantly, the implementation is fast enough for practical applications.

We must recall that our approach is restricted to the correction of point-wise error in structured regions (i.e. base paired nucleotides). Nonetheless it should supplement well existing tools, by using previously discarded information holding, as shown, a strong signal.

Further research, given the potential of error-correction of RNApyro, will evaluate its impact over large datasets with different existing NGS error-correction pipe-line.

Acknowledgments

The authors would like to thank Rob Knight for his suggestions and comments.