Romain Dujardin

SADDLE HYPERBOLICITY IMPLIES HYPERBOLICITY FOR POLYNOMIAL AUTOMORPHISMS OF C 2

Keywords:

We prove that for a polynomial diffeomorphism of C 2 , uniform hyperbolicity on the set of saddle periodic points implies that saddle points are dense in the Julia set. In particular f satisfies Smale's Axiom A on C 2 .

Introduction

Let f be a polynomial automorphism of C 2 with non-trivial dynamics. For such a dynamical system there are two natural definitions for the Julia set. The first one is in terms of normal families: J = J + ∩ J -is the set of points at which which neither (f n ) n≥0 nor (f -n ) n≥0 is locally equicontinuous. The second one is the closure J * of the set of saddle periodic orbits. The inclusion J * ⊂ J is obvious, and whether the reverse inclusion holds is one of the major open questions in higher dimensional holomorphic dynamics.

Following Bedford and Smillie [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF], we say that f is hyperbolic if J is a hyperbolic set for f . Under this assumption we have a rather satisfactory understanding of the global dynamics of f . Indeed it was shown in [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF] that under this assumption the forward and backward Julia sets J + and J -(see §2.1 below for precise definitions) are laminated by stable and unstable manifolds, that the Fatou set is the union of finitely many cycles of attracting basins, that f satisfies Smale's Axiom A on C 2 and finally that J = J * . It was shown by Buzzard and Jenkins [BJ] that f is structurally stable on C 2 . There are also tentative models for a description of the topological dynamics on J (see Ishii [I] for a survey).

On the other hand it is sometimes more natural to postulate that f is uniformly hyperbolic on J * . One reason is that this information can be read off from the periodic points of f . This happens for instance in the study of the stability/bifurcation dichotomy for families of polynomial automorphisms [DL, BD]. The global consequences of hyperbolicity on J * are then less easy to analyze, in particular it does not a priori imply a uniform laminar structure on J ± .

The main result of this paper is that these two notions actually coincide.

Main Theorem. Let f be a polynomial automorphism of C 2 with non-trivial dynamics. If f is hyperbolic on J * , then J = J * .

In particular, if f is hyperbolic on J * , then it is hyperbolic in the sense of [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF].

Recall that the Jacobian Jac(f ) of a polynomial automorphism is a non-zero constant: thus f is dissipative when |Jac(f )| < 1 and conservative when |Jac(f )| = 1. This result was first announced in the dissipative case in [F], but the published proof is not correct 1 , and it has remained an intriguing open problem since then. Recently, Guerini and Peters [GP] managed to establish the result under the more stringent assumption that f is substantially dissipative, that is |Jac(f )| < d -2 , where d is the dynamical degree (see §2.1 for this notion). Observe that only quasi-hyperbolicity on J * is assumed in [GP] while our approach seems to require the full strength of hyperbolicity.

The proof of the main theorem starts with the dissipative case (Section 3). We assume by contradiction that f is dissipative, hyperbolic on J * and that J = J * . In a first stage we show that for some p ∈ J * , J -intersects W s (p) along a non-trivial relatively open subset, which is an unexpected property in the dissipative setting (for instance in the substantially dissipative case, the main point of [GP] is to show that J -∩ W s (p) is totally disconnected). The main input here is the ergodic closing lemma that we obtained in a previous work [START_REF] Dujardin | A closing lemma for polynomial automorphisms of C 2[END_REF].

In a second stage we use the results of [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . VI[END_REF] on the properties of stable slices of J -together with some potential-theoretic ideas to actually derive a contradiction.

The conservative case is treated in Section 4 by a perturbative argument. If f is conservative and hyperbolic on J * , we can find a holomorphic family (f λ ) with f 0 = f containing dissipative parameters, on which J * moves under a holomorphic motion. Again we assume that J * (f ) = J(f ), and use the extension properties of the holomorphic motion of J * obtained in [DL] to derive a contradiction from the previously established dissipative case. Acknowledgment. Thanks to Eric Bedford for many interesting comments.

Preliminaries

In this section we recall some basic facts on the dynamics of polynomial automorphisms of C 2 and hyperbolic dynamics, and establish a few preliminary results.

2.1. Vocabulary and basic facts. Let f be a polynomial diffeomorphism of C 2 with nontrivial dynamics. This is the case exactly when the dynamical degree lim n→∞ (deg(f n )) 1/n is larger than 1. By [FM] there exists a polynomial change of coordinates in which f is expressed as a composition of Hénon mappings (z, w) → (p i (z) + a i w, a i z). We fix such coordinates from now on. The degree of f is d = deg(p i ) ≥ 2 and the relation deg(f n ) = d n holds so that d coincides with the dynamical degree of the original map.

In these adapted coordinates, let

V - R = (z, w) ∈ C 2 , |w| ≥ R, |z| < |w| and V + R = (z, w) ∈ C 2 , |z| ≥ R, |w| < |z| . We fix R 0 > 0 so large that for R ≥ R 0 f (V + R ) ⊂ V + 2R and f -1 (V - R ) ⊂ V - 2R . Hence the points of V + R (resp. V - R ) escape under forward (resp. backward) iteration. We denote by B the bidisk D(0, R 0 ) 2 . The non-wandering set of f is contained in B.
An object (subset, current, or subvariety) in B is said to be vertical (resp. horizontal) if its closure in B is disjoint from {|z| = R 0 } (resp. {|w| = R 0 }). A vertical subvariety has a degree, which is the number of intersection points with a generic horizontal line.

1 A first problem happens in the proof of [START_REF] Fornaess | The Julia set of Hénon maps[END_REF]Thm 2], which corresponds to Step 1 in our proof. Indeed in the construction of the "queer" disk V , the sequence (yn) is contained in W s (J * ) but not a priori in W s loc (J * ), hence one cannot directly deduce that G + (yn) ≥ c. Also, Lemma 6 is not correct: local product structure does not allow to transport whole components of W s (x) ∩ J to components of W s (y) ∩ J when x and y belong to the same global unstable manifold; in particular the boundedness of such a component is not an invariant property.

Here are some standard facts and notation (see e.g. [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . II: Stable manifolds and recurrence[END_REF][START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]):

• K ± is the set of points with bounded forward orbits under f ±1 and K = K + ∩ K -.

Note that K + is vertical in B and f (B ∩ K + ) ⊂ K + . Similarly, K -is horizontal and

f -1 (B ∩ K -) ⊂ K -.
• The complement of K + is denoted by U + and the complement of K -is U -.

• J ± = ∂K ± are the forward and backward Julia sets. If f is dissipative then K -= J -.

• J = J + ∩ J -is the Julia set.

• J * ⊂ J is the closure of the set of saddle periodic points. It is also the support of the unique measure of maximal entropy log d.

The dynamical Green functions G ± are defined by G ± (z, w) = lim d -n log + f ±n (z, w) (where log + (x) = max(log(x), 0)). These are non-negative continuous plurisubharmonic functions on C2 , such that K ± = {G ± = 0} and G ± is pluriharmonic on U +/-:= G +/-> 0 . We let T ± = dd c G ± . The maximum principle implies that Supp(T ± ) = J ± .

The restriction T + | D of T + to a complex submanifold D is a positive measure on D locally defined by ∆(G + | D ) and since G + is continuous this measure coincides 2 with the wedge product

T + ∧ [D]. A useful remark is that if x belongs to K + and D ⊂ C 2 is a holomorphic disk through x along which G + is harmonic, then D ⊂ K + and (f n | ∆ ) n≥1 is a normal family.
If p is a saddle periodic point or more generally if it belongs to a hyperbolic saddle set, it admits stable and unstable manifolds W s/u (p). Each of them is an immersed Riemann surface biholomorphic to C and by [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . II: Stable manifolds and recurrence[END_REF][START_REF] Fornaess | Nessim Complex Hénon mappings in C 2 and Fatou-Bieberbach domains[END_REF]

W s (p) (resp. W u (p)) is dense in J + (resp. J -).
A key point in the present paper is to analyse the topological properties of sets of the form K -∩ W s (p) or K + ∩ W u (p). Following [DL], we define the intrinsic topology to be the topology induced on a stable (resp. unstable) manifold by the biholomorphism W s C, and the corresponding concepts of boundary, interior, etc. will be labelled with the subscript i:

∂ i , Int i , etc.
The following basic lemma will be used several times. By [BLS], the currents T ± have geometric structure, related to the decomposition of J ± into stable and unstable manifolds. By lamination by Riemann surfaces we mean a closed subset L of some open set Ω ⊂ C 2 such that every p ∈ L admits a neighborhood B biholomorphic to a bidisk, such that in the corresponding coordinates, a neighborhood of p in L is a union of disjoint graphs (that is, a holomorphic motion) over the first coordinate in B. A positive current S is uniformly laminar if there is a lamination of Supp(S) by Riemann surfaces and in the corresponding local coordinates S is locally expressed as [∆ a ] dν(a). These disks will be said subordinate to S.

A holomorphic disk D is subordinate to T + if there exists a non-zero uniformly laminar current S ≤ T + such that D is subordinate to S. By [Du1, Prop. 2.3], if p is any saddle point, then any relatively compact disk D ⊂ W s (p) is subordinate to T + . 2.2. Stable (dis)connectivity. It was shown in [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . VI[END_REF] that the connectivity properties of sets of the form K + ∩ W u (p) (resp. K -∩ W s (p)) carry deep information on the geometry of the Julia set. We say that f is stably connected if U + ∩ W s (p) is simply connected for some (and then any) saddle point p, and stably disconnected otherwise. Equivalently, f is stably disconnected if for some saddle point p, W s (p) ∩ K -admits a compact component relative to the intrinsic topology. This actually implies the stronger property that most components of W s (p) ∩ K -are points (see the proof of Lemma 3.2 below for more details).

By [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . VI[END_REF]Cor. 7.4], a dissipative polynomial automorphism is always stably disconnected. It was observed in [START_REF] Dujardin | Some remarks on the connectivity of Julia sets for 2-dimensional diffeomorphisms[END_REF] that this implies a strong non-extremality property for the current T + | B : there exists a decomposition

T + | B = ∞ k=1 T + k where T +
k is an average of integration currents over a family of disjoint vertical disks of degree k (see [START_REF] Dujardin | Some remarks on the connectivity of Julia sets for 2-dimensional diffeomorphisms[END_REF]Thm 2.4]).

Lemma 2.2. Let f be dissipative and hyperbolic on J * and let q ∈ J * . Then q belongs to the support of T + | W u (q) and for (T + | W u (q) )-a.e. q near q, W s B (q ) is a vertical manifold of finite degree in B.

Proof. The first assertion easily follows from the fact that (f n ) n≥0 cannot be a normal family on W u loc (q) (see [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]Lemma 2.8]). The second one is a consequence of [START_REF] Dujardin | Some remarks on the connectivity of Julia sets for 2-dimensional diffeomorphisms[END_REF]Thm 2.4]. Indeed as observed above

T + | B admits a decomposition T + | B = ∞ k=1 T + k where T + k is made of vertical disks of degree k. Thus T + | W u loc (q) = T + ∧ [W u loc (q)] = k T + k ∧ [W u loc (q)]. Now if Γ is a leaf of some T +
k intersecting W u loc (q) at q , then since W u loc (q) is subordinate to T -, q belongs to J * and Γ is a manifold through q along which forward iterates are bounded, hence Γ = W s B (q ). 2.3. Hyperbolicity and local product structure. Let us recall some generalities from hyperbolic dynamics, specialized to our situation. A (saddle) hyperbolic set for f is a compact invariant set Λ ⊂ C 2 such that T C 2 | Λ admits a hyperbolic splitting, i.e. T C 2 | Λ = E s ⊕ E u , where E s and E u are continuous line bundles such that E s (resp. E u ) is uniformly contracted (resp. expanded) by df . Then there exists δ 1 > 0 such that W s δ 1 (Λ) := p∈Λ W s δ 1 (x) and W u δ 1 (Λ) := p∈Λ W u δ 1 (x) form laminations in the δ 1 -neighborhood of Λ. A hyperbolic set is locally maximal if there exists an open neighborhood N of Λ such that Λ = n∈Z f -n (N ). It has local product structure if there exists 0 < δ 2 ≤ δ 1 such that if p, q ∈ Λ are such that d(p, q) < δ 2 then W s δ 1 (p)∩W u δ 1 (q) consists of exactly one point belonging to Λ. It turns out that these two properties are equivalent (see [Y, §4.1]).

We will use the following consequence of the shadowing lemma.

Proposition 2.3. Let Λ be a compact locally maximal hyperbolic set for a polynomial diffeomorphism f of C 2 . Then there exist positive constants η, α and A such that for every n ≥ 0; if x is such that {x, . . . , f n (x)} ⊂ Λ η then there exists y ∈ Λ such that x is Ae -αn -close to the local stable manifold of y.

A similar result holds for negative iterates: if {f -n (x), . . . , x} ⊂ Λ η then there exists z ∈ Λ such that x is Ae -αn -close to the local unstable manifold of z.

The following corollary is well-known.

Corollary 2.4. If Λ is a compact locally maximal hyperbolic set, then

W s (Λ) := x ∈ C 2 , f n (x) -→ n→∞ Λ = p∈Λ W s (p) and similarly W u (Λ) = p∈Λ W u (p).
Note however that W s (Λ), being an increasing union of laminations, doesn't need to have a lamination structure (this is already false when Λ is a hyperbolic fixed point).

Proof of Proposition 2.3 (sketch). This is very classical. Given an orbit segment {x, . . . , f n (x)} as in the statement of the proposition, let y (0) (resp. y (n) ) be a point in Λ such that d(x, y (0) ) < η (resp. d(x, y (n) ) < η). Then define a η-pseudo-orbit (y (k) ) k∈Z as follows

y (k) =      f k y (0) for k < 0; f k (x) for 0 ≤ k ≤ n; f k-n y (n) for k > n.
Then if η is small enough by local maximality and the shadowing lemma there exists a unique y ∈ Λ such that for every k ∈ Z, d f k (y), y (k) < Cη (where C is some constant depending on (f, Λ), see [Y, §4.1]). In particular for 0 ≤ k ≤ n we have d f k (x), f k (y) < Cη and it follows from standard graph transform estimates that d(x, W s δ (y)) ≤ Ae -αn . The next result is a simple application of the techniques of [BLS].

Proposition 2.5. If J * is hyperbolic then it has local product structure. Furthermore global stable and unstable manifolds intersect only in J * :

W s (J * ) ∩ W u (J * ) = J * .
Proof. Hyperbolicity implies that for some δ > 0, if p and q are close enough, W s δ (p) ∩ W u δ (q) consists of a single point r. We have to show that r ∈ J * . Indeed, W s δ (p) (resp. W u δ (q)) is a disk subordinate to T + (resp. T -) so there exists a non-trivial uniformly laminar current S + ≤ T + (resp. S -≤ T -) with W s δ (p) (resp W u δ (q)) as a leaf. By [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]Lem. 8.2], S + and S -have continuous potentials, so the wedge product S + ∧ S -is well defined, and geometric intersection theory of uniformly laminar currents [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]Lem. 8.3] implies that r ∈ Supp(S + ∧ S -). Since S + ∧ S -≤ T + ∧ T -we conclude that r ∈ J * .

The proof of the second assertion is similar. By local product structure,

W s (J * ) = n≥0 f -n (W s δ (J * )),
hence if r ∈ W s (J * ), there exists p ∈ J * such that r ∈ W s (p) so r belongs to a disk subordinate to T + , and likewise r ∈ W u (q) so it belongs to a disk subordinate to T + . Observe that these two disks are distinct: indeed otherwise we would have

W s (p) = W u (q) which is impossible because W s (p) ∩ W u (q) is contained in K which is bounded in C 2 .
So r is an isolated intersection between W s (p) and W u (q) for the leafwise topology. If this intersection is transverse, we argue as above to conclude that r belongs to J * . If it is a tangency, then by [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]Lem. 6.4] for q ∈ J * close to q, we get transverse intersections between W s (p) and W u (q ) close to r and conclude as in the transverse case.

Stability.

A theory of stability and bifurcations for polynomial automorphisms of C 2 was developed in [DL], centered on the notion of weak J * -stability.

A branched holomorphic motion over a complex manifold Λ in C 2 is a family of holomorphic graphs over Λ in Λ × C 2 . It is a holomorphic motion (i.e. an unbranched branched holomorphic motion!) when these graphs are disjoint. A holomorphic family (f λ ) λ∈Λ of polynomial automorphisms of dynamical degree d is weakly J * -stable if the sets J * (f λ ) move under a branched holomorphic motion, and J * -stable if this motion is unbranched. Note that if f λ 0 is uniformly hyperbolic on J * (f λ 0 ), then it is J * -stable near λ 0 in any holomorphic family containing f λ 0 .

A number of properties of weakly J * -stable families are established in [DL], including extension properties of the branched holomorphic motion of J * to K (and more generally to J + ∪ J -), that will be used in Section 4. These properties hold under the standing assumption that the family (f λ ) is substantial 3 : this means that either all members of the family are dissipative, or that no relation of a certain form between multipliers of periodic points persistently holds in the parameter space Λ. Without entering into the details, let us just note that by [START_REF] Buzzard | Kupka-Smale theorem for polynomial automorphisms of C 2 and persistence of heteroclinic intersections[END_REF]Thm 1.4] any open subset of the family of all polynomial automorphisms of dynamical degree d is substantial.

Proof of the main theorem: the dissipative case

The proof is by contradiction so assume that f is a dissipative polynomial diffeomorphism of C 2 , that is uniformly hyperbolic on J * , and that J * J.

Step 1. There exists p ∈ J * and a holomorphic disk

∆ ⊂ W s (p) such that G -| ∆ ≡ 0.
The purpose of the remaining steps 2 and 3 will be to show that such a "queer" component of W s (p) ∩ J -actually does not exist.

Proof. We first claim that it is enough to show that there exists p ∈ J * and q ∈ W s (p) such that q ∈ J \ J * . Indeed observe that W s (p)

∩ J = W s (p) ∩ J -= W s (p) ∩ K -. By Lemma 2.1 we have that ∂ i (W s (p) ∩ J -) ⊂ J * . Hence if q belongs to W s (p) ∩ (J \ J * ), it belongs to the intrinsic interior Int i (W s (p) ∩ J -), hence G -≡ 0 in a neighborhood of q in W s (p).
Let now x ∈ J \ J * . By Corollary 2.4, if ω(x) ⊂ J * then x ∈ W s (J * ) and if α(x) ⊂ J * then x ∈ W u (J * ). Since W u (J * ) ∩ W s (J * ) = J * , we infer that either ω(x) ⊂ J * or α(x) ⊂ J * . In either case we will show that both W s (J * ) ∩ (J \ J * ) and W u (J * ) ∩ (J \ J * ) are non-empty. Thus by symmetry it is enough to deal with the case where ω(x) ⊂ J * .

Choose η so small that Proposition 2.3 holds for J * and ω(x) is not contained in N , where N := (J * ) η is the η-neighborhood of J * .

Consider the sequence of Cesarò averages ν

n = 1 n n-1 k=0 δ f k (x)
. By the ergodic closing lemma of [START_REF] Dujardin | A closing lemma for polynomial automorphisms of C 2[END_REF], every cluster value of the sequence (ν n ) is supported on J * . It follows that the asymptotic proportion of iterates of x belonging to N tends to 1, i.e.

1 n # 0 ≤ k ≤ n -1, f k (x) ∈ N -→ n→∞ 1.
Indeed if a positive proportion of iterates stayed outside N , any cluster limit of ν n would have to give positive mass to N . We thus infer that there are arbitrary long strings {x i , . . . , x i+n } in the orbit of x that are entirely contained in N . Indeed, if on the contrary the length of such a string were uniformly bounded by some n 0 , then the density of iterates outside N would be bounded below by 1/(n 0 + 1). Therefore for every n there exists i n such that {x in , . . . , x in+n } ⊂ N . Choose i n to be minimal with this property. Since ω(x) ⊂ N , there exists j > i such that x j / ∈ N . So finally for every n we can find i n < j n such that j n -i n ≥ n, {x in , . . . , x jn } ⊂ N , x in-1 / ∈ N and x jn+1 / ∈ N .

3 There is an unfortunate terminological conflict here: this should not be confused with the notion of substantial dissipativity mentioned in the introduction.

Let p (resp. p ) be a cluster value of (x in-1 ) (resp. (x jn+1 )). The points p and p belong to J because x does, but not to J * because they lie outside N . It follows from Proposition 2.3 that p ∈ W s (q) for some q ∈ J * and p ∈ W u (q ) for some q ∈ J * . The proof is complete.

Remark 3.1. If f is substantially dissipative i.e. |Jac(f )| < deg(f ) -2
, then the contradiction readily follows from this first step. Indeed Wiman's theorem together with uniform hyperbolicity imply that the vertical degree of components of stable manifolds in some large bidisk B is uniformly bounded (see [START_REF] Guerini | Julia sets of complex Hénon maps[END_REF]Prop. 4.2] or [START_REF] Lyubich | Structure of partially hyperbolic Hénon maps[END_REF]Lem. 5.1]), and it follows that J -∩ W s (x) is totally disconnected for every x (see [START_REF] Dujardin | Some remarks on the connectivity of Julia sets for 2-dimensional diffeomorphisms[END_REF]Thm 2.10] or [START_REF] Guerini | Julia sets of complex Hénon maps[END_REF]Thm. 4.3]), which contradicts the conclusion of Step 1.

In the second and third steps we do not use the assumption that J \ J * = ∅.

Step 2. For every p ∈ J * , if Ω is a component of Int i (W s (p) ∩ J -), then Ω is unbounded for the leafwise topology.

Proof. Note first that by the maximum principle, any component of Int i (W s (p) ∩ J -) = Int i (W s (p) ∩ K -) is simply connected, so Ω is a topological disk. Assume by contradiction that Ω is bounded for the leafwise topology. Then iterating forward a few times if needed, we can suppose that Ω is entirely contained in a local product structure box.

More precisely for small δ > 0, we can fix holomorphic local coordinates (z, w) near p in which p = (0, 0), W s 2δ (p) = {z = 0} and W u 2δ (p) = {w = 0}, and assume Ω is contained in W s δ (p). Note that by Lemma 2.1, ∂ i Ω ⊂ J * . We can assume that for every q ∈ W s δ (p) ∩ J * , W u 2δ (q) contains a graph over the disk D(0, δ) in the first coordinate, with slope bounded by 1/2. Then if |z 0 | ≤ δ, the holonomy h u 0,z 0 along local unstable leaves is well defined on W s δ (p) ∩ J * and maps W s δ (p) ∩ J * into {z = z 0 } ∩ J -. This holonomy is a holomorphic motion so by Slodkowski's theorem [S] it extends to a holomorphic motion of W s δ (p). In particular the motion of ∂ i Ω extends to a motion of Ω and it makes sense to speak about h u 0,z 0 (Ω). This is an open subset of {z = z 0 }, which is topologically a disk and whose boundary is contained in J -. Thus for every n ≥ 0, f -n (∂(h u 0,z 0 (Ω))) in contained in B and by the maximum principle, the same holds for f -n (h u 0,z 0 (Ω)). Finally, O := |z 0 |<δ h u 0,z 0 (Ω) is an open set whose negative iterates remain in B, hence it is contained in the Fatou set of f -1 . But since f is dissipative, this Fatou set is empty, which is the desired contradiction.

Step 3. The unstable holonomy preserves the decomposition

W s (p) = (W s (p) ∩ J -) (W s (p) ∩ U -).
To make this statement precise, observe that for every p ∈ J * , the components of the complement of

∂ i (W s (p) ∩ J -) in W s (p) can be divided into two types: components of Int i (W s (p)∩J -) and components of W s (p)∩U -(note that since U -is open in C 2 , W s (p)∩U -
is open for the intrinsic topology as well). Consider as above local coordinates (z, w) near p in which p = (0, 0), W s 2δ (p) = {z = 0} and W u 2δ (p) = {w = 0}. The unstable holonomy h u 0,z 0 is initially only defined for points of W s δ (p) ∩ J * = ∂ i (W s δ (p) ∩ J -), however by S lodkowski's theorem it can be extended to W s δ (p). By Step 2, components of Int i (W s (p)∩J -) are leafwise unbounded so they cannot be contained in W s 2δ (p). Obviously, the same holds for components of W s (p) ∩ U -.

If q belongs to J * ∩W u δ (p), the extended holonomy h u p,q defines a homeomorphism W s δ (p) → h u p,q (W s δ (p)). By local product structure this homeomorphism preserves J * so any component of W s δ (p) \ J * is mapped onto a component of h u p,q (W s δ (p)) \ J * ), which is itself contained in a component of W s (q) \ J * . The claim of Step 3 is that the extended holonomy h u p,q preserves the type of components.

Since it doesn't make sense to transport a whole leafwise unbounded component by unstable holonomy, to prove this assertion we need to find a criterion that recognizes the type of a component just from local topological properties near a point of its boundary. As already said the maximum principle implies that any component of Int i (W s (p) ∩ J -) is simply connected. Thus Step 3 follows from:

Lemma 3.2. If Ω is a component of W s (p) ∩ U -,
then Ω is not simply connected near any point of ∂Ω, more precisely: if q ∈ ∂Ω and N is any neighborhood of q, there is a loop in N ∩ Ω, homotopic to a point in N , and enclosing a component of W s (p) ∩ J - Proof. Since f is dissipative by [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . VI[END_REF]Cor. 7.4] it is stably disconnected. It follows that almost every unstable component of K + is a point (see [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . VI[END_REF]Thm 7.1] and also [START_REF] Dujardin | Some remarks on the connectivity of Julia sets for 2-dimensional diffeomorphisms[END_REF]Thm 2.10]). More specifically, if µ is the unique measure of maximal entropy, then for µ-a.e. x, the measure T -| W s (x) (which is locally given by the wedge product T -∧ [W s (x)]) gives full mass to the point components of J -∩W s (x). Obviously by Lemma 2.1 every such point component belongs to J * so we can transport it to nearby stable manifolds by unstable holonomy. In addition, the measure T -| W s (x) is holonomy invariant (see [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF]Thm 6.5] or [START_REF] Bedford | Polynomial diffeomorphisms of C 2 . IV. The measure of maximal entropy and laminar currents[END_REF]Thm 4.5

]) so if x is such that T -∧ [W s δ (x)]
gives full mass to point components, then the same holds for nearby x . Thus we conclude that this property holds for every p ∈ J * : T -| W s (p) gives full mass to the point components of J -∩ W s (p).

Lemma 3.3. Let p ∈ J * and Ω be a component of W s (p) ∩ U -such that Ω is locally simply connected near some q ∈ ∂Ω. Then ∂Ω has positive (T -| W s (p) )-measure. This proves Lemma 3.2. Indeed, assuming that Ω is locally simply connected near some q ∈ ∂Ω, Lemma 3.3 asserts that T -| W s (p) carries positive mass on a non-trivial continuum so f cannot be stably disconnected. On the other hand f must be stably disconnected because it is dissipative, and we reach a contradiction.

The idea of Lemma 3.3 is as follows: every neighborhood of q in ∂Ω has positive harmonic measure when viewed from Ω. But the harmonic measure viewed from Ω is absolutely continuous with respect to T -∧ [W s (p)], hence the result. The formalization of this argument requires some elementary potential theory, for which we refer the reader to Doob's classical monograph 4 [Do].

In particular we shall use the formalism of sweeping (or balayage). Let D be a smoothly bounded domain in C, A a non-polar compact subset of D and ν a positive measure on D. The swept measure ρ ν,D,A of ν on A is the distribution on A of the exit point of the Brownian motion in D \ A whose starting point is distributed according to ν. In particular its mass is lower than that of ν since a positive proportion of Brownian paths escape from ∂D. If G ν,D is the Green potential of ν in D, that is the unique negative subharmonic function on D such that G ν,D | ∂D = 0 and ∆G ν,D = ν, then the swept measure of ν on A is ∆R ν,D,A , where (1) R ν,D,A (z) = sup {u(z), u ≤ 0 subharmonic on D and u ≤ G ν,D on A} (see Sections 1.III.4, 1.X and 2.IX.14 in [Do]). If ν and ν have their supports disjoint from A, then the corresponding swept measures are mutually absolutely continuous (as follows from instance from Theorem 1.X.2 in [Do]).

Proof. We first claim that we can shift q slightly so that the assumptions of the lemma hold and in addition W s B (q ) is of bounded vertical degree. Indeed q belongs to J * and for q ∈ W u δ (q) ∩ J * , there is a component of W s δ (q ) \ J * corresponding to Ω under unstable holonomy, which is locally simply connected near q . Since G -is continuous, if q is close enough to q, it takes positive values on that component, so we infer that the property that Ω is a component of U -is open. Now by Lemma 2.2, for (T + | W u δ )-a.e. q , W s B (q ) is a vertical manifold in B of finite degree which establishes our claim. Without loss of generality rename q into q. For every g 0 < min ∂B G -, the component of {G -< g 0 } containing q in W s (q) is relatively compact for the intrinsic topology. We fix such a g 0 which is not a critical value of G -and let D be the corresponding component, which is a smoothly bounded topological disk. From now on we work exclusively in D.

By assumption there is a neighborhood N of q in D and a component U of {G -> 0} ∩ N that is simply connected. We have to show that ∂U ∩N has positive mass relative to ∆G -. We choose N to be closed so that ∂U ∩N is compact. First, observe that for every z 0 ∈ U ∩N , the probability that the Brownian motion issued from z 0 hits ∂U ∩ N before leaving U is positive. Therefore the swept measure ρ δz 0 ,D,∂U ∩N has positive mass and to prove the lemma it is enough to show that is absolutely continuous with respect to ∆G -. Recall that the measure class of the swept measure does not depend on the starting point so we can replace δ z 0 by an arbitrary positive measure on D \ (D ∩ K -). Let 0 < g 1 < g 0 and µ g 1 := ∆(max(G -, g 1 )) be the natural measure induced by G -on the level set {G -= g 1 }. We choose µ g 1 for the initial distribution of Brownian motion. Since G -≡ g 0 on ∂D, the Green function G ∆G -,D of the restriction of ∆G -to D is equal to G --g 0 , and likewise G µg 1 ,D = max(G -, g 1 ) -g 0 .

Thus from (1) we get that

R µg 1 ,D,K -∩D = sup u(z), u s.h. ≤ 0 on D and u ≤ G µg 1 ,D on K -∩ D = sup u(z), u s.h. ≤ 0 on D and u ≤ g 1 -g 0 on K -∩ D = |g 1 -g 0 | G --g 0 g 0 and finally ρ µg 1 ,D,K -∩D = g 0 -g 1 g 0 ∆G -.
The proof is complete.

Step 4. Conclusion.

We just have to assemble the three previous steps. Assume as before by contradiction that f is dissipative, uniformly hyperbolic on J * and J * J. Then by Step 1 there exists p ∈ J * and a "queer" component Ω of W s (p) \ J * along which G -≡ 0. Pick q ∈ ∂Ω. By Lemma 2.1, q ∈ J * so we can follow Ω ∩ W s δ (q) using the holonomy along local unstable manifolds. Then for q ∈ W u loc (q) near q, the holonomy image h q,q (Ω ∩ W s δ (q)) is contained in a queer component of W s (q ) \ J * , which must be leafwise unbounded by Step 2. On the other hand by Lemma 2.2, for generic q in W u loc (q) (relative to the transverse measure T + | W u loc (q) ) W s B (q ) is of bounded degree, in particular any component of K -∩ W s B (q ) is leafwise bounded. This contradiction finishes the proof. 4. Proof of the main theorem: the conservative case Again the proof is by contradiction, so assume that f is a conservative polynomial automorphism of C 2 such that J * is a hyperbolic set and J * J. We will use a perturbative argument and the dissipative case of the theorem to reach a contradiction.

Assume that f is written as a product of Hénon mappings f = h 1 • • • • • h k and let (f λ ) λ∈B be a parameterization of a neighborhood of f in the space of such products, that is, the space of coefficients of the h i , and such that f 0 = f . We can assume that B is a ball in C N for some N . Since λ → Jac(f λ ) is an open map, there exist parameters arbitrary close to 0 for which f λ is dissipative. As already said, by [START_REF] Buzzard | Kupka-Smale theorem for polynomial automorphisms of C 2 and persistence of heteroclinic intersections[END_REF]Thm 1.4] there is no persistent relation between multipliers of periodic orbits so the family is substantial in the sense of [DL].

Since f 0 is hyperbolic on J * (f 0 ) the family (f λ ) is J * -stable in a neighborhood of the origin, that is, J * (f λ ) moves under a holomorphic motion. Reducing the parameter space we can assume that J * is hyperbolic throughout B. Pick a point p = p(0) ∈ J(f 0 ) \ J * (f 0 ). It was shown in [START_REF] Dujardin | Stability and bifurcations for dissipative polynomial automorphisms of C 2[END_REF]Thm 5.12] that in a (weakly) J * -stable family, the motion of J * extends to a branched holomorphic motion of K. Thus there exists a holomorphic continuation p(λ) of p(0) such that for every λ ∈ B, p(λ) belongs to K(f λ ). Furthermore for every λ ∈ B, p(λ) is disjoint from J * (f λ ). Indeed if for some λ 0 ∈ B we had p(λ 0 ) ∈ J * (f λ 0 ), then by [START_REF] Dujardin | Stability and bifurcations for dissipative polynomial automorphisms of C 2[END_REF]Lem 4.10] p(λ) would have to coincide throughout the family (f λ ) with the natural continuation of p(λ 0 ) as a point of the hyperbolic set J * , which is not the case since p(0) / ∈ J * (f 0 ). Let now λ 1 ∈ B be such that f λ 1 is dissipative. Then by the first part of the proof J(f λ 1 ) = J * (f λ 1 ), and K(f λ 1 ) \ J(f λ 1 ) is non-empty since it contains p(λ 1 ). For a dissipative hyperbolic map K \ J = (K + ∩ J -) \ (J + ∩ J -) = Int(K + ) ∩ J -, so we deduce that Int(K + (f λ 1 )) is non-empty. By [START_REF] Bedford | Polynomial diffeomorphisms of C 2 : currents, equilibrium measure and hyperbolicity[END_REF], Int(K + (f λ 1 )) is a finite union of attracting basins of periodic sinks, therefore f λ 1 admits an attracting periodic point. On the other hand by [START_REF] Dujardin | Stability and bifurcations for dissipative polynomial automorphisms of C 2[END_REF]Thm 4.2], periodic points stay of constant type in a J * -stable family (this holds even in the presence of conservative maps, provided the family is substantial), so f 0 must have an attracting orbit, which is contradictory since it is conservative. The proof is complete.

  Lemma 2.1 ([DL, Lemma 5.1]). Let p be a saddle periodic point. Then the boundary of W s (p) ∩ J -relative to the intrinsic topology in W s (p) is contained in J * . We denote by W s B (p) the connected component of W s (p) ∩ B containing p (and accordingly for W u ). Likewise, W s δ (p) is the connected component of W s (p) ∩ B(p, δ) containing p, and W s loc (p) denotes an unspecified open neighborhood of p in W s (p).

It is standard to define d c = i 2π (∂ -∂).Accordingly, ∆ here is 1/2π times the ordinary Laplacian.

Research partially supported by ANR project LAMBDA, ANR-13-BS01-0002 and a grant from the Institut Universitaire de France.


PhD Fanny Orlhac
email: orlhacf@gmail.com

MD, PhD Augustin Lecler

MD Julien Savatovski

PhD Jessica Goya-Outi

PhD Christophe Nioche

MD Frédérique Charbonneau

PhD Nicholas Ayache

PhD Frédérique Frouin

MD Loïc Duron

PhD Irène Buvat

How can we combat multicenter variability in MR radiomics? Validation of a correction procedure

Keywords: Magnetic Resonance Imaging, Neoplasms / diagnostic imaging, Image Processing, Computer-Assisted / methods* CT: computed tomography CE-T1w: contrast enhanced T1-weighted D1/D2: prostate cancer patient database 1/2 GANs: generative adversarial networks GLCM: gray-level co-occurrence matrix hWS: hybrid white stripe LDA: Linear Discriminant Analysis MRI: magnetic resonance imaging PET: positron emission tomography

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

1. Radiomic feature values obtained using different MR scanners or imaging protocols can be harmonized by combining off-the-shelf image standardization and feature realignment procedures. 2. Harmonized radiomic features enable one to pool data from different scanners and centers without substantial loss of statistical power caused by intra-and inter-center variability. 3. The proposed realignment method is applicable to radiomic features from different MR sequences and tumor types and does not rely on any phantom acquisition.

Introduction

Radiomics faces the critical issue of lack of reproducibility that still hampers the successful translation of radiomic model discovery into better diagnosis, patient classification or monitoring radiomicsbased tools. Indeed, radiomic features values are significantly affected by the technical settings of the imaging devices and protocols, as demonstrated in Positron Emission Tomography (PET), Computed Tomography (CT), and Magnetic Resonance Imaging (MRI) [START_REF] Yan | Impact of image reconstruction settings on texture features in 18F-FDG PET[END_REF][START_REF] Berenguer | Radiomics of CT features may be nonreproducible and redundant: influence of CT acquisition parameters[END_REF][START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF]. To tackle the variability of radiomic features induced by different technical settings, radiomic models can be designed using a wide variety of images encompassing most technical settings, or image or radiomic feature values have to be harmonized before designing models. Ignoring the center effect, as is often observed in many papers, results in lack of generalization of the radiomic models [START_REF] Reuzé | Prediction of cervical cancer recurrence using textural features extracted from 18F-FDG PET images acquired with different scanners[END_REF].

In prospective studies, imaging protocols could be harmonized upstream between centers to minimize the impact of imaging protocols on feature values [START_REF] Boellaard | FDG PET/CT: EANM procedure guidelines for tumour imaging: version 2.0[END_REF][START_REF] Clarke | The Quantitative Imaging Network: NCI's historical perspective and planned goals[END_REF], although harmonizing between machines of different generations often comes with degrading the image quality achieved by the most recent scanners [START_REF] Boellaard | FDG PET/CT: EANM procedure guidelines for tumour imaging: version 2.0[END_REF]. In retrospective studies, this approach is not an option. Several groups have proposed to reduce the variability by resampling the images to a common voxel size or by filtering the images to match spatial resolution [START_REF] Shafiq-Ul-Hassan | Voxel size and gray level normalization of CT radiomic features in lung cancer[END_REF][START_REF] Mackin | Harmonizing the pixel size in retrospective computed tomography radiomics studies[END_REF]. However, this requires accessing the images retrospectively and the filtering procedure reduces the quality of images acquired using the most recent devices. Others apply a z-score transformation [START_REF] Chatterjee | Creating robust predictive radiomic models for data from independent institutions using normalization[END_REF] to each feature value based on mean and standard deviations measured in each center for that feature, but this assumes that images produced by the different centers have been obtained in similar patient samples (eg, same proportion of advanced and early stage tumors), which is sometimes difficult to achieve.

In genomics, researchers face a similar problem called batch effect and caused by the handling of samples by different laboratories, different technicians, on different days that can obscure individual variations. To deal with that problem in genomics, Johnson et al [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF] introduced the ComBat realignment method. The method realigns all data in a single space in which the batch effect is discarded without altering the biological information. This approach has already been successfully validated for radiomic features measured from PET [START_REF] Orlhac | A postreconstruction harmonization method for multicenter radiomic studies in PET[END_REF] and CT [START_REF] Orlhac | Validation of a method to compensate multicenter effects affecting CT radiomics[END_REF][START_REF] Mahon | ComBat harmonization for radiomic features in independent phantom and lung cancer patient computed tomography datasets[END_REF] images of patient or phantom data in studies supporting the relevance of harmonization.

In MR, the challenge is even more difficult as, unlike in PET and CT where images are expressed in kBq/mL and Hounsfield Units respectively, there is no standard MR intensity grey scale, implying the lack of a tissue-specific absolute intensity numeric meaning, even within the same MR imaging protocol, for the same body region, for images obtained on the same scanner, for the same patient.

The standardization of image intensities among patients is therefore absolutely needed for comparing values of intensity-based features. In brain MR, standardization approaches have been proposed to correct for the intensity variability [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF][START_REF] Zhuge | Intensity standardization simplifies brain MR image segmentation[END_REF][START_REF] Ge | Numerical tissue characterization in MS via standardization of the MR image intensity scale[END_REF][START_REF] Nyúl | On standardizing the MR image intensity scale[END_REF]. In particular, the hybrid white stripe (hWS) method proved to be successful in the context of neurodegenerative diseases and cancer [START_REF] Shinohara | Statistical normalization techniques for magnetic resonance imaging[END_REF][START_REF] Kickingereder | Radiogenomics of glioblastoma: machine learning-based classification of molecular characteristics by using multiparametric and multiregional MR imaging features[END_REF]. ComBat has been validated in MRI for the harmonization of cortical thickness measurements across scanners [START_REF] Fortin | Harmonization of cortical thickness measurements across scanners and sites[END_REF].

Although it has been used in MR radiomic studies [START_REF] Lucia | External validation of a combined PET and MRI radiomics model for prediction of recurrence in cervical cancer patients treated with chemoradiotherapy[END_REF][START_REF] Whitney | Harmonization of radiomic features of breast lesions across international DCE-MRI datasets[END_REF][START_REF] Wang | Preoperative MRI-based radiomic machine-learning nomogram may accurately distinguish between benign and malignant soft-tissue lesions: a two-center study[END_REF][START_REF] Zhang | Pretreatment MRI radiomics analysis allows for reliable prediction of local recurrence in non-metastatic T4 nasopharyngeal carcinoma[END_REF], it has never been validated in that highly challenging context.

Here, we extend this approach to provide a harmonization procedure applicable to any radiomic feature. We demonstrate that by combining the image standardization (such as hWS) with ComBat realignment, MR radiomic features can be pooled without being adversely impacted by multiple sources of variability, ensuring higher sensitivity and specificity of multicenter MR radiomic studies.

Materials and Methods

The study was carried out in accordance with the World Medical Association's Declaration of Helsinki.

For experiment 2 involving MRI brain studies, the institutional review board of the Fondation Ophtalmologique A. Rothschild approved the study (IRB No. 1512-016-726), and the requirement to obtain written informed consent was waived because of the retrospective nature of the study. For experiment 3 using MRI prostate studies, all patient data are publicly available in https://datadryad.org/resource/doi:10.5061/dryad.b3d257g and were initially studied in [START_REF] Penzias | Identifying the morphologic basis for radiomic features in distinguishing different Gleason grades of prostate cancer on MRI: preliminary findings[END_REF]. All patient data were anonymized. All authors had control of the data and information submitted for publication.

Experiment 1: phantom studies

The phantom data used in this study have been extracted from the RIDER Phantom MRI study [START_REF] Jackson | Magnetic resonance assessment of response to therapy: tumor change measurement, truth data and error sources[END_REF] and are publicly available in the TCIA platform [START_REF] Clark | The Cancer Imaging Archive (TCIA): maintaining and operating a public information repository[END_REF]. The phantom consists of 19 doped gel filled tubes containing a gadolinium-based contrast agent. We used the T1 acquisitions obtained with a 1.5T scanner and a 3T scanner (respectively scanners B and D in [START_REF] Jackson | Magnetic resonance assessment of response to therapy: tumor change measurement, truth data and error sources[END_REF], details in Supplemental data 1). For each image, 19 spherical Volumes of Interest (VOI) of 3.5 mL centered on each tube were drawn. We computed 42 radiomic features (Supplemental data 2) using LIFEx freeware [START_REF] Nioche | LIFEx: a freeware for radiomic feature calculation in multimodality imaging to accelerate advances in the characterization of tumor heterogeneity[END_REF] (www.lifexsoft.org), including an open-source radiomic protocol compliant with the Image Biomarker Standardisation

Initiative guidelines [START_REF] Zwanenburg | The Image Biomarker Standardization Initiative: standardized quantitative radiomics for high throughput image-based phenotyping[END_REF]. Radiomic features were calculated using a fixed bin size [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF][START_REF] Orlhac | 18F-FDG PET-derived textural indices reflect tissue-specific uptake pattern in non-small cell lung cancer[END_REF] set to the average standard-deviation of the signal intensity, between the minimum and the maximum intensity measured in all VOI. This discretization step is required to set voxels with similar intensity to the same value hence to reduce the impact of noise.

Experiment 2: MRI brain studies

For experiment 2, we retrospectively selected 18 patients (13 men; mean age, 50±18 years; age range, 26-85 years; Table 1) with grade III and IV glial tumors from January 2017 to May 2018 from an institutional database. All patients underwent two MRI scans using the same protocol: one on a at 1.5T scanner (Philips Achieva, Philips Medical Systems) and the other one on a 3T scanner (Philips Ingenia). The median delay between the two scans was 30 days (range: [4-93 days]) without chemotherapy, surgery, radiotherapy and any visual evolution of the tumor and tumor heterogeneity between the scans. Two MR sequences (details in Supplemental data 1) were acquired: a 3D FLAIR ( 17patients) and a 3D contrast enhanced T1-weighted (CE-T1w) scan (14 patients).

For each patient and each sequence, the 3T images were coregistered to the 1.5T images using rigid transformations in FSL-FLIRT [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]. Field inhomogeneity was corrected using the N4 algorithm [START_REF] Tustison | N4ITK: improved N3 bias correction[END_REF] owing to the publicly available ANTs software (http://stnava.github.io/ANTs) with the standard setting of hyper-parameters.

For each sequence, the tumor lesions were manually segmented based on a consensus of two radiologists (A.L. and L.D. with 9 and 2 years of experience, respectively) on the 1.5T images and the resulting regions were copied on the 3T images. Three slices (top, middle, bottom) were selected in each tumor to obtain three 2D-regions of interest (ROI) per tumor, yielding a total of 54 tumor ROI for FLAIR images (=3x18 tumors; one patient had two distinct lesions) and 51 tumor ROI for CE-T1w images (=3x17 tumors; one patient had two distinct lesions and another had three). In addition, in each patient, 6 regions of 0.5 mL each were drawn in the white matter (WM), yielding 102 WM VOI for FLAIR images and 84 WM-VOI for CE-T1w images that were copied onto the 3T images.

Each patient image volume was standardized irrespective of the other patients using the hWS method [START_REF] Shinohara | Statistical normalization techniques for magnetic resonance imaging[END_REF] as previously described [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF]. The hWS method applies a z-score transformation to the brain voxel values based on the normal-appearing WM intensities distribution.

For each ROI and VOI based on native and hWS-standardized images resampled at 1x1x1 mm 3 , we computed 42 radiomic features using LIFEx. Radiomic features were calculated using a fixed bin size [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF][START_REF] Orlhac | 18F-FDG PET-derived textural indices reflect tissue-specific uptake pattern in non-small cell lung cancer[END_REF] set to the average standard-deviation of the WM signal intensity, between the minimum and the maximum intensity measured in all WM and tumor VOI for each sequence separately (details in Supplemental data 1).

Experiment 3: MRI prostate studies

Two prostate cancer patient databases (D1 and D2; Table 1) with publicly available radiomic features were used [START_REF] Penzias | Identifying the morphologic basis for radiomic features in distinguishing different Gleason grades of prostate cancer on MRI: preliminary findings[END_REF]. These two databases have been initially entirely studied by an independent team to investigate the relationship between features computed from MR images and from digitized tissue images in order to discriminate between prostate cancer grades, without taking into account that MRI scans were acquired in two different centers [START_REF] Penzias | Identifying the morphologic basis for radiomic features in distinguishing different Gleason grades of prostate cancer on MRI: preliminary findings[END_REF]. Here, we precisely investigate how accounting for the center effect actually changes the ability of each MRI feature to distinguish between tumor grades. 1). Based on a co-registration with histology images, the corresponding tumor regions were manually segmented by a radiologist. MR images were standardized to a template distribution based on the per-patient median of intra-prostatic pixel intensities of D1 [START_REF] Nyúl | New variants of a method of MRI scale standardization[END_REF]. For each region, 2379 radiomic features were computed using a home-made software (details of feature calculation described in [START_REF] Penzias | Identifying the morphologic basis for radiomic features in distinguishing different Gleason grades of prostate cancer on MRI: preliminary findings[END_REF]) and we selected the 2326 features available for all patients for our analysis.

Realignment method

To correct for the scanner effect, the ComBat realignment method was used [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF]. In the context of radiomics, ComBat has already been validated for PET [START_REF] Orlhac | A postreconstruction harmonization method for multicenter radiomic studies in PET[END_REF] and CT features [START_REF] Orlhac | Validation of a method to compensate multicenter effects affecting CT radiomics[END_REF][START_REF] Mahon | ComBat harmonization for radiomic features in independent phantom and lung cancer patient computed tomography datasets[END_REF]. The method directly applies to the radiomic feature values and estimates the scanner-effect by matching the statistical distributions of the feature values measured in VOI j for each scanner i:

y ij =α+𝑋 𝑖𝑗 β+ γ i +δ i ε ij
where α is the average value for feature y ij , 𝑋 is a design matrix for the covariates of interest; 𝛽 is the vector of regression coefficients corresponding to each covariate, γ i is an additive scanner effect, and δ i is a multiplicative protocol effect affected by an error term (ε ij ). The model parameters α, 𝛽, γ i and δ i are estimated using a maximum likelihood approach based on the set of available observations from the two scanners in experiments 1 and 2 and based on the two patient databases for experiment 3.

The corrected values are obtained using:

y ij ComBat = y ij -α ̂-𝑋 𝑖𝑗 𝛽 ̂-γ i δi ̂+α ̂
where α ̂, 𝛽 ̂, γ i ̂ and δ i ̂ are estimators of α, 𝛽, γ i and δ i .

The non-parametric form of the model was used, with no assumption regarding the statistical laws followed by the features and a transformation determined for each feature separately. For experiments 1 and 2, no biologic covariate was used (i.e. 𝑋=0) since the data came from the same patients or phantom scanned on 1.5T and 3T machines, and we realigned feature values computed from WM and tumor regions in patient data separately. For experiment 3, we introduced the Gleason grade as a binary covariate since the proportion of low vs intermediate/high-risk regions was very different between the 2 databases (32% low-risk VOI in D1, 68% in D2, Table 1).

To facilitate the access to the ComBat method for medical imaging professionals, we provide a free on-line application (available at https://forlhac.shinyapps.io/Shiny_ComBat/), named ComBaTool, with example input files (Supplemental data 3-4) and a step-by-step tutorial (Supplemental data 5). This application embeds a free function called ComBat [START_REF] Fortin | Harmonization of cortical thickness measurements across scanners and sites[END_REF] (https://github.com/Jfortin1/ComBatHarmonization) based on the R software but running the application does neither require R or any third-party software to be installed, nor to have any programming skills.

Statistical analysis

Statistical analysis was performed with the R software (version 3.6.1).

In experiment 1, we performed univariate two-sided Friedman tests before and after ComBat realignment between the two phantom scans. In experiment 2, we used two-sided Friedman tests for each radiomic feature to test whether the values derived from the 1.5T and 3T scans were significantly different both in the WM and in the tumor regions in three configurations: C1) native images without ComBat realignment; C2) hWS-standardized images without realignment; C3) hWS-standardized images with realignment. The Benjamini-Hochberg procedure was used to control the false discovery rate [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]. P values less than 0.05 were interpreted as statistically significant. Bland-Altman graphs were plotted to demonstrate the differences in feature values calculated from the 1.5T and 3T scans.

In

Results

Patient characteristics are shown in Table 1.

Experiment 1

In ), demonstrating that the scanner effect was no longer detectable for the vast majority of radiomic features. Figure 1 shows the evolution of the distribution of the Correlation radiomic feature calculated from the gray-level co-occurrence matrix (GLCM) . On native FLAIR images, the plot shows a shift in distribution with greater values for WM-VOI and tumor lesions for 3T scans compared to 1.5T scans. After hWS standardization and realignment, the distributions between the two scanners better overlap. To clarify the respective role of hWS and ComBat, Figure 2 shows the Bland-Altman plots of the mean value measured in WM-VOI for FLAIR images based on 3T scans and 1.5T scans. The hWS standardization within each patient rescaled the values to make them similar between the two scans. The realignment reduced the systematic difference between the two.

The same trends were observed for CE-T1w images (Table 2; Supplemental data 7).

Experiment 3

On T2w prostate images after standardization performed by [START_REF] Penzias | Identifying the morphologic basis for radiomic features in distinguishing different Gleason grades of prostate cancer on MRI: preliminary findings[END_REF], 461 out of 2326 radiomic features had P values of Wilcoxon tests less than 0.05 for distinguishing between low and intermediate/high risk when pooling the two patient cohorts (D1+D2). After ComBat without any co-variate, 460 out of 2326 P values were less than 0.05. Using the Gleason grade co-variate in ComBat, 636 out of 2326 P values were less than 0.05. Figure 3 demonstrates a better alignment of radiomic feature values extracted from low-risk VOI and intermediate/high-risk VOI separately between the two patient groups after using ComBat with a co-variate accounting for the recruitment specificity of each center.

When a risk (low or moderate/high) was randomly assigned to each VOI, no P value was less than 0.05 before and after ComBat without and with a co-variate representing the Gleason grade.

The multivariate radiomic model identified using LDA on the D1 data to distinguish low versus intermediate/high risk was applied to D2 patients, yielding a Youden Index of 0.12 (Sensitivity=19%, Specificity=93%) before ComBat. After ComBat, the Youden Index increased to 0.20 (Sensitivity=27%, Specificity=93%) and to 0.43 (Sensitivity=58%, Specificity=86%) using the Gleason grade as co-variate in ComBat.

Discussion

The scanner effect affects the radiomic feature values extracted from MR images, introducing major confounding factors in multicentric or multi-protocol studies. Here, we validated a harmonization procedure combining ComBat realignment with MR-image standardization to co-analyze MR radiomic features extracted from different scanners. Using phantom data and brain scans acquired for the same patients (without any tumor evolution detected visually between the two scans) with 1.5T and 3T scanners, we showed that this harmonization procedure realigns radiomic feature distributions and removes the scanner effect for T1, FLAIR and CE-T1w images. The goal was not to test our ability to reproduce feature values measured in 3T MR images from 1.5T images, since we expect different signals from the two devices with more details in the 3T images (cf Figure 1). Yet, in the context of radiomics, pooling images acquired using different devices, different acquisition and reconstruction protocols is often needed to increase the size of cohort. In that context, we demonstrated that ComBat could realign feature values so that all data could be analyzed together, even if images had been acquired with different magnetic fields. It is important to underline that a different ComBat transformation is estimated for each sequence and each tissue type independently because imaging protocols do not have the same effect on each tissue. Using the prostate scans acquired in different patients from two centers, we confirmed the effectiveness of the harmonization for T2w images and demonstrated that harmonization did not alter the discriminant information conveyed by the features. This experiment also shows that pooling data corrected for the scanner effect could increase the statistical power, identify more radiomic features able to distinguish between the low-risk and intermediate/high-risk regions in prostate lesions and yield a more discriminant multivariate model.

Importantly, we showed that when no difference between groups was expected, here between the sham low-risk and intermediate/high-risk VOI, ComBat did not introduce any false positive differences.

The ComBat realignment method is fast, easy to use and operates directly on radiomic feature values (no training set needed, no phantom acquisition, no need to access images). It is applicable to radiomic features extracted from different MR sequences after a first step of image standardization, as previously described [START_REF] Goya-Outi | Computation of reliable textural indices from multimodal brain MRI: suggestions based on a study of patients with diffuse intrinsic pontine glioma[END_REF]. We also demonstrated the added value of the covariate in the realignment process when patient characteristics are different between centers (here Gleason grade) for univariate and multivariate analyses. To deal with the center effect, other authors reported the potential of Generative Adversarial Networks (GANs) to transform images from one imaging protocol (or a domain) to another [START_REF] Qu | Wavelet-based semi-supervised adversarial learning for synthesizing realistic 7T from 3T MRI[END_REF]. Although promising results have been reported in the literature [START_REF] Zhong | Inter-site harmonization based on dual generative adversarial networks for diffusion tensor imaging: application to neonatal white matter development[END_REF][START_REF] Modanwal | MRI image harmonization using cycleconsistent generative adversarial network[END_REF], these techniques require access to the images, unlike ComBat. The ComBat realignment method has been previously used in MR radiomic studies [START_REF] Lucia | External validation of a combined PET and MRI radiomics model for prediction of recurrence in cervical cancer patients treated with chemoradiotherapy[END_REF][START_REF] Whitney | Harmonization of radiomic features of breast lesions across international DCE-MRI datasets[END_REF][START_REF] Wang | Preoperative MRI-based radiomic machine-learning nomogram may accurately distinguish between benign and malignant soft-tissue lesions: a two-center study[END_REF][START_REF] Zhang | Pretreatment MRI radiomics analysis allows for reliable prediction of local recurrence in non-metastatic T4 nasopharyngeal carcinoma[END_REF] without any explicit validation or investigation of the respective role of the image standardization and of the scanner/protocol effect compensation as studied here (Figures 1 and2). In [START_REF] Lucia | External validation of a combined PET and MRI radiomics model for prediction of recurrence in cervical cancer patients treated with chemoradiotherapy[END_REF], authors reported an increased accuracy of Entropy extracted from apparent diffusion coefficient MR images to predict the locoregional control in cervical cancer after ComBat, fully consistent with our findings.

Our study has some limitations. We could only include 18 patients in Experiment 2 because it is very uncommon for patients to undergo MR both on 1.5 and 3T scanners within a time lapse during which the tumor has not visually evolved. Still, this small sample allowed us to confirm results obtained using the phantom data. In addition, such a small number allowed us to demonstrate that ComBat performed well even with a limited number of cases, confirming results published in genomic applications [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF]. Another limitation is that our findings should still be validated for other cancer types, MR sequences and devices.

In conclusion, we demonstrated that the ComBat realignment method in combination with intrapatient image standardization could efficiently remove the scanner/protocol effect while preserving the individual variations in phantom, brain and prostate MR scans. This approach enables large MR multicentric studies to investigate the added value of radiomic analysis in patient management. To facilitate large multicenter/multi-protocol radiomic studies, we provide the ComBat method as an online ComBaTool application.

Figure 1 :

 1 Figure 1: Experiment 2: 18 patients with brain lesions were scanned on both 1.5T and 3T scanners.Based on native or for hybrid White Stripe (hWS)-standardized, 42 radiomic features were computed in a tumor region and in a white matter region. As an example, the probability density function (%) of the Correlation radiomic feature calculated from the gray-level co-occurrence matrix (GLCM) on FLAIR images is plotted here without and with ComBat realignment (ComBaTool was applied separately on the two tissue types: white matter and tumor) for 1.5T MRI (in orange) and 3T MRI (in blue). P values are for Friedman tests of each tissue between the two MRI devices.

Figure 2 :

 2 Figure 2: Experiment 2: Bland-Altman plots of the mean value computed in white matter regions based on 1.5T and 3T scans for FLAIR native images (A), for hybrid White Stripe (hWS)-standardized images (B) and for hWS-standardized images with ComBat realignment (C).

Figure 3 :

 3 Figure 3: Experiment 3: boxplots of Feature #20 (called Gabor:cos:theta=0:lambda=2:Standard Deviation in [24]) for low-risk VOI and intermediate/high-risk VOI, before ComBat realignment (A, D), after ComBat realignment without covariate (B, E) and after ComBat realignment with covariate (C, F) for the prostate patient cohorts D1 and D2 separately (A, B, C) or together (D, E, F). P values are from Wilcoxon tests.

  experiment 3, we performed Wilcoxon tests with the Benjamini-Hochberg procedure for all

	radiomic features to distinguish between low and intermediate/high-risk groups when pooling
	patients from D1 and D2, without ComBat realignment, with realignment and with realignment
	including the Gleason grade as a covariate of interest. To show that ComBat does not create false
	positive results, we repeated these tests after randomly assigning a label to each VOI to get 53 sham
	low-risk VOI and 52 sham intermediate/high-risk VOI. To identify the risk group, we built a multivariate
	signature by means of a Linear Discriminant Analysis (LDA) using D1 dataset as a training set and
	including only the features with a p-value of univariate Wilcoxon test less than 5%. We tested the
	classification performance on D2 data by calculating the Youden Index (=Sensitivity+Specificity-1). We
	repeated this procedure in three configurations: without ComBat realignment, with realignment and
	with realignment including the Gleason grade as a covariate of interest.

Table 1 :

 1 patient characteristics. Mean age of patients is not reported in [Penzias et al. PlosOne 2018]. CE-T1w = Contrast Enhanced T1-weighted; D1 = prostate dataset 1; D2= prostate dataset 2; T2w = T2-weighted.

	Parameter	MRI brain studies	MRI prostate studies
	Sex
	Men	13	D1 = 23; D2 = 13
	Women	5	0
	Mean age (y)	50±18	Not reported
	No. of tumor regions
	FLAIR images	18 3D-regions (54 2D-regions)	-
	CE-T1w images	17 3D-regions (51 2D-regions)	-
	T2w images	-	D1 = 65; D2 = 40
	Gleason grade
	Low risk	-	D1 = 21(32%); D2 = 26 (65%)
	Intermediate/high risk	-	D1 = 44 (68%); D2 = 14 (35%)
	Note. --

Table 2 :

 2 significant Friedman tests (P<0.05) without and with hybrid White Stripe (hWS) standardization and/or ComBat realignment for brain scans.

		w/o hWS	with hWS	with hWS
		standardization	standardization	standardization
		w/o ComBat	w/o ComBat	with ComBat
	FLAIR images
	WM	37/42 (88)	29/42 (69)	0/42 (0)
	Tumor lesions	41/42 (98)	25/42 (60)	1/42 (2)
	CE-T1w images
	WM	37/42 (88)	27/42 (64)	0/42 (0)
	Tumor lesions	36/42 (86)	2/42 (5)	0/42 (0)

Note. --Data are numerator/denominator; data in parentheses are percentage. WM = White Matter. CE-T1w = Contrast Enhanced T1-weighted. w/o = without.

Supplemental data legends

Supplemental data 1: MR acquisition and radiomic analysis parameters for experiments 1 and 2.

Supplemental


Julien Fontchastagner
email: julien.fontchastagner@ensem.inpl-nancy.fr

Frédéric Messine
email: frederic.messine@enseeiht.fr

Yvan Lefèvre
email: yvan.lefevre@enseeiht.fr

An Exact Optimization Code Combined with a Hybrid Model for Magnetic Couplings Design

The purpose of this paper is to show a new methodology for solving inverse problems of design. It is based on the use of a deterministic global optimization algorithm. In front of limitations induced by the use of analytical models, the authors present a way to associate Interval Branch and Bounds techniques and numerical computations based on finite element method. Thanks to this methodology, we are able to solve exactly the associated inverse problem of a magnetic coupling design.

I. INTRODUCTION

Magnetic couplings are very useful devices for many applications as seal-less pumps in the chemical and petrochemical industries, or the aeronautical and maritime ones for instance. Their functionalities are to transmit a motion between two separated zones without mechanical contacts. The studied structure is a co-axial synchronous magnetic torque coupling. It consists of two rings of permanent magnets separated by an insulating partition. On each side of it, we find an airgap and a binding band.

A. Previous Methodology

In [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Optimization Code[END_REF], the authors solve the problem of such devices design by applying a rational methodology. The problem is understood and defined as an inverse problem, i.e. from some characteristic values given by the schedule of conditions (for example the torque), get the dimensions as well as the structure and the composition of a co-axial magnetic torque coupling. These inverse problems are formulated as a mixed constrained global optimization problem.

In order to solve exactly the so-formulated global optimization problems, a particular algorithm (called IBBA) based on a Branch and Bound technique where the bounds are computed using interval analysis has been developed and extended. See [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Optimization Code[END_REF], [START_REF] Messine | Essays and Surveys in Global Optimization, Chap.10 "A Deterministic Global Optimization Algorithm for Design Problems[END_REF] and [START_REF] Fontchastagner | Design of Electrical Rotating Machines by Associating Deterministic Global Optimization Algorithm With Combinatorial Analytical and Numerical Models[END_REF] for details.

The advantages of these kinds of methods are :

• they use mixed variables (real, integer, boolean, . . . );

• they do not need a starting point (or a set of ones);

• they guarantee to obtain the global minimum of the problem (they are deterministic and global).

Their main drawbacks are their computation times and the fact that they need (until nowadays) explicit analytic expressions of criteria and constraints.

The inverse problem of couplings design problem is formulated as a mixed constrained global optimization problem defined as follows :

           min x∈I R nr , z∈I N ne , b∈B n b , σ∈ nc i=1 K i f (x, z, σ, b) g i (x, z, σ, b) ≤ 0 ∀i ∈ {1, . . . , n g } Γ(x, z, σ, b) = Γ f ixed (1)
where f is a real function (for example the volume), K i represents an enumerated set of categorical variables (for instance, the kind of magnet), and B = {0, 1} is the boolean set (for the fact that there are binding bands or not). g i are some geometrical constraints. Moreover, an equality constraint upon the maximum electromagnetic torque Γ is added, traducing the fact that it must be equal to a fixed value (Γ f ixed ). The torque expression comes from an analytical model based on the resolution of Poisson's equations using the separation of variables by keeping only the term due to the first non-zero harmonic (the fundamental) [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Optimization Code[END_REF].

Several techniques of constraints propagation, and others adaptations were applied in order to reduce the time convergence. Indeed, it was very difficult to compute efficient bounds with interval arithmetic, because of the complexity of our torque formula. Thanks to these techniques, the authors obtained encouraging results dealing with the minimization of the global volume or with the volume of magnets. In our knowledge, it was the first time that such a problem was solved.

B. Need of a More Accurate Model

If a finite element software is used to compute the electromagnetic torque and the mean value of flux density in yokes (the physical constraints of our problem), we notice a nonnegligible difference between the values given by our model and the numerical ones (around 20%). However, the kind of analytical models used is known to provide results close to those obtained by a numerical method [START_REF] Couderc | Analitycal modelling of high speed electrical machines in view of their optimal design[END_REF]. Actually, the observed gap is due to the fact that we do not have considered the harmonics greater than the fundamental, and moreover the optimization process is inclined to maximize the error due to these assumptions.

II. NEW APPROACH : USE A HYBRID MODEL

The idea is to associate the advantages of two kinds of models : the swiftness of analytics and the accuracy of finite elements based techniques.

The analytic expression of the electromagnetic torque is given by the resolution of Poisson's equations using the separation of variables. Instead of taking into account only the fundamental (Γ o1 ) as in [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Optimization Code[END_REF], the used expression (Γ o3 ) is the sum of the two first non-zero harmonics (1 and 3). The thicknesses of iron yokes are deduced from the Gauss's law of magnetism and the fact that the mean value of flux density in the yokes B y must be less or equal to the maximum value B M (σ y ) above which the iron is definitively saturated.

A specially dedicated finite element code has been written in order to automatically draw the geometry, create the mesh, and perform a magneto-static resolution of a magnetic coupling. In output, the torque Γ FE (computed with the Maxwell Stress Tensor) and the mean values of flux density in the inner and outer yokes (B FE yi and B FE yo ) are given in less than 1 second. This tool is in the form of a black-box (named NUMTFD) which can be easily called by another program.

III. FORMALIZATION OF THE NEW DESIGN PROBLEM

Now the question is to know how to associate such a model with an Interval Branch and Bound Algorithm.

A. Inverse Problem Formulation

The design problem must be re-formulated. In fact it is not possible to directly include our black-box as a constraint. Indeed, when the optimization process begins, the lengths of the intervals are too large to perform a valid finite element computation [START_REF] Fontchastagner | Design of Electrical Rotating Machines by Associating Deterministic Global Optimization Algorithm With Combinatorial Analytical and Numerical Models[END_REF].

The idea is to define several zones corresponding to different expressions. If we are far from the wanted value (Γ f ixed ), we use Γ o1 . If we get close, but not close enough to call Γ FE , Γ o3 is used. The associate problem is then :

                                 min x∈I R nr , z∈I N ne , b∈B n b , σ∈ nc i=1 K i f (x, z, σ, b) g j (x, z, σ, b) ≤ 0 ∀j ∈ {1, . . . , n g }      (1 -λ ana )Γ f ixed ≤ Γ o1 ≤ (1 + λ ana )Γ f ixed (1 -λ FE )Γ f ixed ≤ Γ o3 ≤ (1 + λ FE )Γ f ixed Γ FE (x, z, σ, b) = Γ f ixed B FE yi (x, z, σ, b) ≤ B M (σ yi ) B FE yo (x, z, σ, b) ≤ B M (σ yo )
where λ ana and λ FE are used to define the different zones (for the tests, we have chosen 40% and 20%). Two inequality constraints are added to ensure that the yokes are not saturated.

B. New Interval Branch & Bound Algorithm

The corresponding Interval Branch & Bound Algorithm (IBBA+NUMTFD) has been coded. Its principle is to bisect the initial domain into smaller and smaller boxes and then to eliminate the boxes where the global optimum cannot occurs:

• by proving, using interval bounds, that no point in a box can produce a better solution than the current best one; • by proving, using interval arithmetic, that at least one constraint cannot be satisfied by any point in such a box. Some techniques of constraint propagation and limitations have been included to improve the convergence as in [START_REF] Fontchastagner | Some Co-axial Magnetic Couplings Designed Using an Analytical Model and an Exact Global Optimization Code[END_REF].

IV. RESULTS AND CONCLUSION

The problem corresponding to the minimization of the global volume of a studied magnetic coupling is solved and results are given in Table I (strike-through text represents nonsatisfied constraints). We have chosen Γ f ixed = 10 N•m, the real parameters represent geometric quantities, the integer p is the number of poles pairs, and the different σ values are used to choose the kinds of materials. To compare with the previous methodology, the problem is solved using only Γ o1 (IBBA), next using Γ o3 (IBBAo3), and finally using our hybrid model (IBBA+NUMTFD).

Only the results corresponding to the new methodology answer perfectly to our non-homogeneous mixed constrained global optimization problem.


Isabelle Ydens

Philippe Degée

Cécile Nouvel

Edith Dellacherie

Jean-Luc Six

Philippe Dubois
email: philippe.dubois@umh.ac.be

'Surfactant-free' stable nanoparticles from biodegradable and amphiphilic poly(ε-caprolactone)-grafted dextran copolymers

A variety of poly(ε-caprolactone)-grafted dextran copolymers were synthesized with a controlled architecture through a three-step procedure: partial protection of the dextran hydroxyl groups by silylation, ring-opening polymerization of ε-caprolactone initiated from remaining free hydroxyl groups on partially silylated dextran, and silyl ether deprotection under mild conditions. The potential of these amphiphilic grafted copolymers as surfactants was investigated by water/toluene interfacial tension measurements. Since each copolymer showed noticeable surfactant properties, their ability to form 'surfactant-free' nanoparticles was evaluated using the nanoprecipitation method. It was found that remarkably stable core (polyester) -shell (polysaccharide) nanoparticles with a mean diameter around 200 nm were formed.

Introduction

The development of drug delivery systems aimed at optimizing drug release and targeted therapy has contributed to significant advancements in pharmaceutical engineering of novel formulations such as nanoparticles, which are solid colloidal polymeric carriers less than 1 µm in size. Several publications highlighted the ability of nanoparticles to reduce adverse effects of various drugs released at non-target sites, typically hydrophobic ones [1][2][3][4][5]. Most commonly used methods for preparing polymer-based nanoparticles include emulsion-evaporation [6], in situ monomer polymerization [7], a method based on the salting out effect [8][9], and nanoprecipitation [10,11]. The latter originally developed by Fessi et al. [11] represents an easy and reproducible technique, which has been quickly and widely explored by several research groups for producing both vesicle and matrix type nanoparticles, also termed nanocapsule and nanosphere, respectively [10,[12][13]. The main disadvantages of conventionally formed nanoparticles are the need for large amounts of surfactants such as poly(vinyl alcohol) [14,15], Tween ® , etc. [16,17], to stabilize the suspension, and the difficulty to remove these tensioactive agents at the end of the production process. Moreover, a lot of properties including biodegradability, size distribution and time-evolution as well as drug release depend on the concentration of residual surfactants [18][19][20].

For the last decade, surfactant-free nanoparticles formation has been investigated by several authors [11,[21][22][23]. Fessi et al. developed surfactant-free nanocapsules of poly(D,L-lactide) (PLA) by the nanoprecipitation technique [11]. More recently, poly-(lactide-co-glycolide) nanoparticles were prepared by using poly(ethylene glycol)based block copolymers as substitutes for conventional surfactants [23]. Dialysis post-treatment was also used for producing small and narrowly distributed nanoparticles using conventional surfactants or amphiphilic block and graft copolymers [24][25][26]. The present paper aims at reporting on the exclusive use of amphiphilic, biodegradable poly(ε-caprolactone)-grafted dextran (Dex-g-PCL) enriched in hydrophobic poly(ε-caprolactone) (PCL) for preparing surfactant-free nanoparticles through the nanoprecipitation technique. Purposely, several well-defined Dex-g-PCL compositions were synthesized according to a three-step procedure already reported by us and depicted in Scheme 1 [27,28]. Beside their macromolecular characterization, the interfacial properties of these amphiphilic graft copolymers were studied by the pendant drop method, as well as the mean diameter and size distribution of the soformed nanoparticles by using dynamic light scattering.

Partial and reversible silylation

Scheme 1. Three-step synthetic pathway to dextran-g-poly(ε-caprolactone) (Dex-g-PCL) copolymers

Experimental part

Materials

Dextrans were purchased from Pharmacia Biotech (M n = 6600 or 21 300; M w /M n = 1.8). Dex-g-PCLs were synthesized according to an original and controlled three-step procedure recently reported by us and illustrated in Scheme 1 [27,28]. Macromolecular characteristics of the Dex-g-PCL graft copolymers used in the present study are shown in Tab. 1 using the following acronyms: Dex x -gy PCL z , where x and z denote the molar masses of dextran and poly(ε-caprolactone) grafts (in kg/mol), respectively, and y denotes the average number of PCL grafts per 100 glucose units.

Dextran-g-poly(ε-caprolactone) copolymer (Dex x -gy PCL z ) synthesis [27,28] The first step consisted of the partial reversible silylation of dextran hydroxyl groups in order to control and limit the number of free hydroxyl groups and to facilitate polysaccharide solubilization in low polar organic solvents suitable to perform the controlled ring-opening polymerization (ROP) of ε-caprolactone (CL). We showed that high silylation yields could be reached without any degradation of dextran by using hexamethyldisilazane in dimethyl sulfoxide (DMSO) while hydroxyl groups located on the third and fourth carbon atoms of glucose repeating units were by far less reactive [30][31]. The second step thus consisted in the regioselective initiation of CL ROP from the hydroxyl groups remaining free along the highly silylated dextran backbone, and previously activated by the addition of triethylaluminium in catalytic amount. The high efficiency of grafting and the control over the PCL grafts' molar mass and molar mass distribution was evidenced by NMR and size exclusion chromatography (SEC), and attested the 'living' character of the coordinationinsertion mechanism of the CL ROP. The final step involved the selective hydrolytic deactivation of the propagating active species, together with the removal of the trimethylsilyl protective groups under mild acidic conditions, which did not affect dextran and PCL molecular parameters.

Nanoparticle preparation

Dex-g-PCL graft copolymers were dissolved in DMSO (99+% from Aldrich) at 50°C until complete solubilization (1 h). The organic solutions with concentrations ranging from 0.5 to 2 mg/mL were then added dropwise with vigorous magnetic stirring into a defined volume excess of Millipore water previously thermostated at 50°C (0.03 ≤ DMSO/H 2 O (v/v) ≤ 0.11). The resulting opalescent colloidal suspension was then dialyzed against Millipore water for 5 days using a 3500 g/mol molecular cut-off dialysis tubular membrane (Spectra/Por® Biotech RC membranes). Finally, the particles suspension was concentrated by water volatilization to a final volume of 9 mL.

Characterization

Interfacial tensions were determined according to the pendant drop method using a Drop Shape Analysis System DSA 10-MK2 tensiometer (from Wilten Fysika) equipped with a thermostated chamber and a circulator Thermo Haake DC10. Practically, a drop of Millipore water was formed at the capillary tip and immersed into a quartz cell containing various Dex-g-PCL copolymer solutions in toluene (99+%, from Devos-François) with concentrations ranging from 10 -5 to 10 g/L. Stock amphiphilic graft copolymer solutions were prepared at room temperature by dissolving a known amount of Dex-g-PCL in the appropriate volume of toluene (conc. = 10 g/L). Solutions of lower concentration were obtained by subsequent dilution of the stock solutions. Interfacial tension was determined at 20°C after a time long enough to reach constant values from the shape of the pendant drop by fitting the Gauss-Laplace equation to the experimental drop shape coordinates [29]. All measurements were triplicate and displayed a maximum variation lower than 2%. Dynamic light scattering measurements were carried out at 25°C using a BI-160 apparatus (Brookhaven Instruments Corporation, USA) with He-Ne laser source operating at 17 mW and delivering a vertically polarized light (λ = 633 nm). The particle size and size distribution were calculated using the CONTIN algorithm. For each sample, the mean diameter and standard deviation were calculated over three measurements.

Tab. 1 H NMR spectroscopy [27]. c Number-average molar mass of the PCL grafts as determined by gravimetry and 1 H NMR spectroscopy [27]. z denotes the molar mass in kg/mol. d PCL weight fraction in the copolymer as determined by gravimetry and confirmed by 1 H NMR spectroscopy [27].

Results and discussion

Macromolecular characteristics of studied dextran-g-poly(ε-caprolactone) copolymers

Tab. 1 shows the main macromolecular characteristics of the investigated dextran-gpoly(ε-caprolactone) graft copolymers, Dex x -gy PCL z , where x and z denote dextran and poly(ε-caprolactone) (PCL) number-average molar masses (in kg/mol), respectively, and y denotes the average number of PCL grafts per 100 glucose units. These amphiphilic and biodegradable graft copolymers were synthesized according to a recently reported three-step procedure illustrated in Scheme 1 [27][28]30]. The solubility of the so-formed Dex x -gy PCL z graft copolymers was shown to strongly depend on the relative content in PCL grafts, i.e., the weight fraction in PCL (F wPCL ). For instance, limited solubility in water was observed when F wPCL was higher than c. 0.3. PCL-enriched graft copolymers were rather soluble in organic solvents such as toluene or tetrahydrofuran. When solubilized in DMSO at c. 50°C, which is a common solvent for both native dextran and PCL segments, Dex x -gy PCL z graft copolymers with F wPCL > 0.3 were suspected to adopt a core-shell conformation with an inner stretched hydrophilic polysaccharide chain surrounded by an outer shell of hydrophobic polyester grafts. 1 H NMR spectroscopy in d 6 -DMSO gave credit to such a conformation with a screening of polysaccharide protons resulting from a restriction of dextran chain mobility and a decrease of relaxation times [27]. In this further study, it has been decided to investigate the amphiphilic behaviour of these graft copolymers and their ability to generate stable nanoparticles in aqueous medium.

Amphiphilic character of Dex x -gy PCL z graft copolymers

In order to get a better insight into the potential of Dex x -gy PCL z graft copolymers to behave as surfactants, interfacial tensions were measured by using the pendant drop method applied to drops of water dipped into toluene copolymer solutions of various concentrations. The semi-logarithmic plot of the interfacial tension (γ in mN/m) vs. concentration (in g/L) was typically sigmoidal whatever the studied copolymer composition (Fig. 1). Critical concentrations (C crit ) were determined from the intersection between the tangents drawn from higher concentration portions of the sigmoidal plots [32]. For the different copolymers studied, C crit is ranging from 0.3 to 2.3•10 -3 g/L (Tab. 2). Critical concentration tends to increase with the content in hydrophobic component (F wPCL ) at least for a given dextran molar mass (entries 1 and 2, Tab. 2). It also comes out that increasing dextran molar mass, while keeping PCL weight fraction and average number of polyester grafts (y) almost constant, contributes to lower C crit (entries 1 and 3, Tab. 2). Finally, higher molecular area (A) and effectiveness (Π crit = γ crit -γ 0 [33]) are reached for Dex 6,6 -g-76 PCL 0,2 , i.e., for the copolymer with the highest average number of PCL grafts per 100 glucosidic units, which is in good agreement with observations from other similar amphiphilic systems [34,35]. Interestingly enough, a similar behaviour has been recently reported for dextran-gpolylactide copolymers (Dex-g-PLA) [37]. It was observed that these copolymers showed noticeable surfactant properties whatever their polyester weight fraction. Moreover, depending on their solubility, the self-organization of Dex-g-PLA copolymer in water or toluene was evaluated using fluorescence spectroscopy. It was proved that these amphiphilic grafted copolymers were able to produce either hydrophobic or hydrated microdomains in water or toluene solutions, respectively. As a conclusion, these studies confirm that the PCL-grafted dextran copolymers can behave as surfactants at the interface between two immiscible solvents such as toluene and water and prompt us to investigate their ability to form nanoparticles by nanoprecipitation. c Critical concentration determined from the intersection between the tangents drawn from higher concentration portions of the sigmoidal plots of γ vs. log C/(g•L -1 ). d Interfacial tension at the critical concentration.

Tab. 2. Effect of PCL-grafted dextran copolymer composition on the main interfacial parameters

Entry Dex x -g-y PCL z a F wPCL b C crit c /(mg•L -1 ) γ crit d /(mN•m -1 ) A e /nm 2 Π crit f /(mN•m -1 )
e Molecular area calculated from the following equations:

C T R . log d Γ 303 2 γ d ex ⋅ ⋅ ⋅ ⋅ - = ; N A ⋅ = ex Γ 1
where R is the gas constant, T the absolute temperature, Γ ex is the surface excess of surfactant molecules (in mol/area), C is the surfactant concentration, and N is Avogadro's number. f Π crit = γ crit -γ 0 , where γ 0 and γ crit stand for the interfacial tensions of the water drop immersed in pure toluene and a copolymer solution in toluene at critical concentration, respectively.

Nanoparticle formation

As sketched in Fig. 2, Dex x -gy PCL z copolymers were first dissolved in DMSO at 50°C until complete solubilization. The organic solution was then added dropwise and with vigorous magnetic stirring into a large volume excess of water previously thermostated at 50°C (see Exptl. part). The suspension rapidly turned opalescent due to colloidal nanoparticle formation. From the aforementioned surface tension properties it can be assumed that the water-soluble dextran backbone is preferably orientated towards the aqueous phase, whereas the hydrophobic PCL side chains concentrate in the dispersed organic phase. To separate DMSO from the medium, the suspension was dialyzed against distilled water using a 3500 g/mol molecular cut-off dialysis tubular membrane. Finally, the resulting suspension was concentrated to a final volume of 9 mL. The influence of the initial DMSO-to-H 2 O volume ratio (R) and of the initial copolymer concentration in DMSO on the size distribution of nanoparticles was investigated by dynamic light scattering at 25°C. It is worth mentioning that no accurate determination of solution viscosity and differential refractive index was performed so that only apparent mean diameters of nanoparticles can be approached. Solutions of selected dextran-g-PCL graft copolymers were prepared in DMSO (1 mg/mL) and defined volumes were added into 9 mL of Millipore water at 50°C with vigorous magnetic stirring. As shown in Fig. 3, increasing R from 0.03 to 0.11 contributed to increase the apparent mean diameter of the so-generated nanoparticles. Except the specific case of the nanoparticles prepared from Dex 6,6 -g-52 PCL 1,9 at high DMSO-to-water volume ratio (R = 0.11) that display a bimodal size distribution, all other studied compositions yielded monodisperse nanoparticles. In other words, these experiments demonstrate that the smaller the DMSO relative volume compared to water the better defined the nanoparticles' distribution. From Fig. 4, the reproducibility of the nanoprecipitation technique can be appreciated since two individually prepared batches of Dex 6,6 -g-76 PCL 0,2 copolymer-based nanoparticles exhibit similar mean diameter and size distribution for R = 0.11.

Effect of initial copolymer concentration in DMSO and stability of colloidal suspension

This series of experiments was exclusively realized on Dex 6,6 -g-76 PCL 0,2 copolymer as it displayed a monomodal size distribution whatever the R value studied, at least in our prevailing experimental conditions. Practically, 1 mL of Dex 6,6 -g-76 PCL 0,2 copolymer solutions in DMSO with concentrations ranging from 0.5 to 2 mg/mL was added to 9 mL of pure water. It comes out that no significant effect of copolymer concentration on the mean diameter can be detected (Fig. 5). Furthermore, the time dependence of the mean size of the recovered nanoparticles has been investigated as a tool to evaluate their stability. Remarkably enough, nanoparticles are stable at least within 3 weeks when stored at room temperature. Indeed, an apparent mean diameter of 201 ± 19 nm is measured after three weeks for Dex 6,6 -g-76 PCL 0,2 based nanoparticles (R = 0.11), which is close to the starting value of 188 ± 16 nm (Fig. 6). It can be concluded that dextran backbones likely located at the surface of nanoparticles are able to efficiently stabilize PCL-enriched nanoparticle suspensions in pure water. Interestingly, these observations are in agreement with recently reported data by Gref et al. [36] on other PCL (low content) graft copolymer based nanoparticles prepared by the 'interfacial migration/solvent evaporation' method. There is no doubt that such behaviour paves the way to the formulation of drug nanocarriers for controlled release and modulated biodistribution.

Conclusion

Well-defined PCL-grafted dextran amphiphilic copolymers have shown high effectiveness in reducing water/toluene interfacial tension. Extending the length of both dextran backbone and PCL grafts allows reducing the critical concentration for a given weight composition. In a next step, these dextran-g-PCL copolymers have been used to readily prepare 'surfactant-free' nanoparticles by the nanoprecipitation technique in water. Whatever the composition of PCL-enriched graft copolymers, adequate experimental conditions in terms of copolymer concentration in DMSO and DMSO-to-H 2 O volume ratio can be found allowing formation of nanoparticles with an apparent mean diameter close to 200 nm and narrow size distribution. As a rule, the higher the PCL content, the lower should be the relative DMSO volume compared to pure water while the initial concentration of the graft copolymer in DMSO showed no significant effect, at least for concentrations ranging from 0.2 to 0.5 mg/mL. From the time dependence of the nanoparticles' mean diameter as recorded over a period of three weeks, it can be concluded that stable nanoparticle suspensions have been obtained. Such behaviour is consistent with the formation of a core-shell structure consisting of a hydrophobic polyester core surrounded by a hydrophilic dextran shell, which prevents rapid agglomeration.

Fig. 1 .

 1 Fig. 1. Semi-logarithmic plot of the concentration dependence of the interfacial tension (γ) for various dextran-g-poly(ε-caprolactone) compositions (Tab. 1)

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Illustration of the nanoprecipitation technique applied to dextran-g-poly(εcaprolactone) copolymers

5 Fig. 5 .

 55 Fig. 5. Influence of Dex 6,6 -g-76 PCL 0,2 copolymer initial concentration in DMSO on the apparent mean diameter of nanoparticles as determined by dynamic light scattering (DMSO-to-H 2 O volume ratio R = 0.11)

Fig. 6 .

 6 Fig. 6. Apparent size distribution of nanoparticles obtained from Dex 6,6 -g-76 PCL 0,2 as determined by light scattering (initial copolymer concentration in DMSO 1 mg/mL, and DMSO-to-H 2 O volume ratio R = 0.11)

  Number-average molar mass of dextran as determined by SEC in water at 35°C with reference to dextran standards. x denotes the molar mass in kg/mol.

	1. Macromolecular characteristics of dextran-g-PCL amphiphilic graft copoly-
	mers (Dex x -g-y PCL z )
	Entry	Sample	M n Dex /(g•mol -1 ) a	y b M n PCL /(g•mol -1 ) c	F wPCL	d
	1	Dex 6,6 -g-52 PCL 1,9	6600	52	1900	0.87
	2	Dex 6,6 -g-76 PCL 0,2	6600	76	200	0.52
	3	Dex 21,3 -g-48 PCL 2,9	21 300 48	2900	0.89

a b Average number of PCL grafts per 100 glucosidic units as determined by

  Graft copolymer acronyms, where x and z denote the molar masses of dextran and PCL grafts in kg/mol, respectively, and y denotes the average number of PCL grafts per 100 glucosidic units. b PCL weight fraction in the copolymer as determined by gravimetry and confirmed by 1 H NMR spectroscopy[27].

	1	Dex 6,6 -g-52 PCL 1,9	0.87	2.3	27.3	0.5	9.3
	2	Dex 6,6 -g-76 PCL 0,2	0.52	1.3	27.1	0.9	11.0
	3	Dex 21,3 -g-48 PCL 2,9	0.89	0.3	27.6	0.5	8.3
	a

Acknowledgement: This work was partially supported by both the Région Wallonne and Fonds Social Européen in the frame of Objectif 1-Hainaut: Materia Nova program. I. Y. is grateful to F.R.I.A. for her Ph.D. grant. This work was also partially supported by a 'Tournesol' Program 2000 -2001 (Program N° 00827WH). LPCM thanks the Belgian Federal Government Office of Science Policy (SSTC-PAI 5/3) for general support.