Anestis Antoniadis
email: anestis.antoniadis@imag.fr

Xavier Brossat -Jairo

Cugliari -Jean-Michel Poggi

Mai

Xavier Brossat
email: xavier.brossat@edf.fr

Jairo Cugliari
email: jairo.cugliari@math.u-psud.fr

Jean-Michel Poggi
email: jean-michel.poggi@math.u-psud.fr

Optimization, Learning and Statistical Methods

Keywords: functional data, clustering, electricity power demand, wavelets

We present two methods for detecting patterns and clusters in high dimensional time-dependent functional data. Our methods are based on wavelet-based similarity measures, since wavelets are well suited for identifying highly discriminant local time and scale features. The multiresolution aspect of the wavelet transform provides a time-scale decomposition of the signals allowing to visualize and to cluster the functional data into homogeneous groups. For each input function, through its empirical orthogonal wavelet transform the first method uses the distribution of energy across scales generate a handy number of features that can be sufficient to still make the signals well distinguishable. Our new similarity measure combined with an efficient feature selection technique in the wavelet domain is then used within more or less classical clustering algorithms to effectively differentiate among high dimensional populations. The second method uses dissimilarity measures between the whole time-scale representations and are based on wavelet-coherence tools. The clustering is then performed using a k-centroid algorithm starting from these dissimilarities. Practical performance of these methods that jointly designs both the feature selection in the wavelet domain and the classification distance is demonstrated through simulations as well as daily profiles of the French electricity power demand.

Classification non supervise des donnes fonctionnelles l'aide d'ondelettes

Résumé : Nous prsentons deux mthodes de dtection des clusters dans des donnes fonctionnelles avec une structure temporelle. Nos mthodes sont bases sur des mesures de similarit fonds dans la dcomposition en ondelettes, car cette dcomposition est bien adapte pour identifier des caractristiques localises en temps et chelles. L'aspect multi-rsolution de la transforme en ondelettes permet une dcomposition des signaux en temps et chelle permettant de visualiser et de regrouper les donnes fonctionnelles en groupes homognes. La premire mthode permet d'obtenir, travers la transform orthogonal empirique en ondelettes, la distribution de l'nergie travers les chelles et ainsi de gnrer un faible nombre d'attributs encore capables de bien distinguer les signaux. La nouvelle mesure de similarit associe une technique efficace de slection des variables dans le domaine des ondelettes est ensuite utilise dans des algorithmes de clustering classiques. La deuxime mthode utilise des mesures de dissimilitude entre les reprsentations temps-chelle bass sur la notion de cohrence en ondelettes. Nous utilisons un algorithme de k-centroides pour obtenir les groupes partir de ces dissimilitudes. La performance pratique de nos mthodes est dmontre travers de simulations et des profiles de consommation journalire d'nergie lectrique en France.

Mots-clés : donnes fonctionnelles, classification non supervise, demande d'nergie lectrique, ondelettes

Introduction

In different fields of applications, explanatory variables are not standard multivariate observations, but are functions observed either discretely or continuously. Ramsay and Dalzell (1991) gave the name "functional data analysis" to the analysis of data of this kind. As evidenced in the work by Ramsay andSilverman (1997, 2002) (see also Ferraty and Vieu (2006)), a growing interest is notable in investigating the dependence relationships between complex functional data such as curves, spectra, time series or more generally signals. Functional data often arise from measurements on fine time grids, and if the sampling grid is sufficiently dense, the resulting data may be viewed as a sample of curves. These curves may vary in shape, both in amplitude and phase. Typical examples involving functional data can be found when studying the forecasting of electricity consumption, temporal gene expression analysis or ozone concentration in environmental studies to cite only a few.

Given a sample of curves, an important task is to search for homogeneous subgroups of curves using clustering and classification. Clustering is one of the most frequently used data mining techniques, which is an unsupervised learning process for partitioning a data set into sub-groups so that the instances within a group are similar to each other and are very dissimilar to the instances of other groups. In a functional context clustering helps to identify representative curve patterns and individuals who are very likely involved in the same or similar processes. Recently, several functional clustering methods have been developed such as variants of the k-means method (Tarpey and Kinateder (2003); Tarpey (2007); Cuesta-Albertos and Fraiman (2007)) and clustering after transformation and smoothing (CATS) (Serban and Wasserman ( 2004)) to model-based procedures, such as clustering sparsely sampled functional data (James and Sugar ( 2003)) or mixed effects modeling approach using B-splines (Luan and Li (2003)) that mostly concentrate on curves exhibiting a regular behaviour.

Our interest in time series of curves is motivated by an application in forecasting a functional time series when the most recent curve is observed. This situation arises frequently when a seasonal univariate continuous time series is segmented into consecutive segments, for example days, and treated as a discrete time series of functions. The idea of forming a functional valued discrete time series from segmentation of a seasonal univariate time series has been introduced by Bosq (1991). Suppose one observes a square integrable continuous time stochastic process X = (X(t), t ∈ R) over the interval [0, T ], T > 0 at a relatively high sampling frequency. The commonly used approach is to divide the interval [0, T ] into sub-intervals [(l -1)δ, lδ], l = 1, . . . , n with δ = T /n, and to consider the functional-valued discrete time stochastic process Z = (Z i , i ∈ N), defined by

Z i (t) = X(t + (i -1)δ), i ∈ N ∀t ∈ [0, δ). (1)
The random functions Z i thus obtained, while exhibiting a possibly nonstationary behavior within each continuous time subinterval, form a functional discrete times series that is usually assumed to be stationary. Such a procedure allows to handle seasonal variation of size δ in a natural way. This set-up has been used for prediction when ones consider a Hilbert-valued discrete time stationary autoregressive processes (see Bosq (1991); Besse and Cardot (1996); Pumo (1992); Antoniadis and Sapatinas (2003)) or for more general continuous-time processes (see Antoniadis et al. (2006)). However, as already noticed above, for many functional data the segmentation into subintervals of length δ may not suffice to make reasonable the stationary hypothesis of the resulting segments, that is the key for efficient prediction. For instance, in modeling the electrical power demand process the seasonal effect of temperature and the calendar configuration strongly affects the mean level and the shape of the daily load demand profile. Recognizing this, our aim is therefore to propose a clustering technique that clusters the functional valued (discrete) times series segments into groups that may be considered as stationary so that in each group more or less standard functional prediction procedures such as the one cited above can be applied.

We will apply the methodology focusing on edf's ( Électricité de France1 ) national power demand for a year. This is essentially a continuous process even though we only count with discrete records sampled at 30 minutes for the whole year. Some of the facts associated to the electricity power demand induce to think that the process is not stationary. We will construct a functional data set by splitting the continuous process as in equation ( 1) where the parameter δ will be a day.

Although slicing an univariate time series produces functional data, we do not observe the whole segments but a sample of the values at some time points. One could then use a vector representation of each observation. For example Wang et al. (2008) proposed to measure the distance between observations through the high dimensional multivariate distribution of all sampled time points along each curve. This approach does not exploit the eventually potential information of correlations between points of a single curve. To avoid this, many authors have clustered the coefficients of a suitable basis representation of functions. Since the analyzed curves are infinite-dimensional and temporalstructured, one projects each curve over an appropriate basis of the functional space to extract specific features form the data which are then used as inputs for clustering or classification. One may cite for example Abraham et al. (2003) where the authors proposed to cluster the spline fit coefficients of the curves using k-means, or James and Sugar (2003) that use a spline decomposition specially adapted for sparsely sample functional data. Nevertheless, attention must be paid to the chosen basis because this operation involves linear transformation of data that may not be invariant for the clustering technique (see Tarpey ( 2007)). Splines are often used to describe functions with a certain degree of regularity. However, we will be working with curves like in Figure 7 that may present quite irregular paths. We chose to work with wavelets because of their good approximations properties on sample paths that might be rather irregular. One may note here that similar methodologies relying upon wavelet decompositions for clustering or classifying time series have been developed in the literature (e.g. see Pittner et al. (1999) and Ray & Mallick (2006)).

Wavelets offer an excellent framework when data are not stationary. For example, in Gurley et al. (2003) the wavelet transform is used to develop the concept of wavelet-coherence that describes the local correlation of the timescale representation of two functions. Grinsted et al. (2004) proved that this concept is convenient for clustering geophysical time series. Another example supporting such a fact is the work by Quiroga et al. (2004) which uses wavelets to detect and cluster spikes on neural activity. Motivated by this and the fact that the wavelet transform has the property of time-frequency localization of the time series, we propose hereafter a time-series feature extraction algorithm using orthogonal wavelets for automatically choosing feature dimensionality for clustering. We also study some more complex alternatives using the waveletcoherence concept in sake of better exploiting the well localized information of the wavelet transform.

The rest of the paper is organized as follows. Section 2 is a reminder on multiresolution analysis and introduces the basis supporting our feature extraction algorithm by means of the energy operator. Following wavelet analysis we cluster the functional data using the extracted features in Section 3. Our first clustering algorithm uses k-means as unsupervised learning routine. We test the proposed method in Section 4 on simulated and real data. Section 5 presents a more sophisticated method for clustering functional data using a more specific dissimilarity measure. Finally, we conclude the paper by summarizing the main contributions and perspectives in Section 6.

Feature extraction with wavelets

In this section we first introduce some basic ideas of the wavelet analysis before introducing more specific material: the energy contributions of the scale levels of the wavelet transform which are the key tools for future clustering. More details about wavelets and wavelet transforms can be found for example in Mallat (1999).

We will consider a probability space (Ω, F , P ) where we define a functionvalued random variable Z : Ω → H, where H is a (real) separable Hilbert space (e.g. H = L 2 (T ) the space of squared-integrable functions on T = [0, 1) (finite energy signals) or H = W s 2 (T ) the Sobolev space of s-smooth function on T , with integer regularity index s ≥ 1) endowed with the Hilbert inner product < ., . > H and the Hilbert norm . H .

Wavelet transform

A wavelet transform (WT for short) is a domain transform technique for hierarchical decomposing finite energy signals. It allows a real valued function to be described in terms of an approximation of the original function, plus a set of details that range from coarse to fine. The property of wavelets is that the broad trend of the input function is captured in the approximation part, while the localized changes are kept in the detail components. For short, a wavelet is a smooth and quickly vanishing oscillating function with good localization properties in both frequency and time. This is suitable for approximating curves that contain localized structures. A compactly supported WT uses a orthonormal basis of waveforms derived from scaling (i.e. dilating or stretching) and translating a compactly supported scaling function φ and a compactly supported mother wavelet ψ. We consider periodized wavelets in order to work over the interval [0, 1], denoting by

φ(t) = l∈Z φ(t -l) and ψ(t) = l∈Z ψ(t -l), for t ∈ [0, 1], RR n°7515
the periodized scaling function and wavelet, that we dilate or stretch and translate φ j,k (t) = 2 j/2 φ(2 j t -k), ψ j,k (t) = 2 j/2 ψ(2 j t -k).

For any j 0 ≥ 0, the collection

{φ j0,k , k = 0, 1, . . . , 2 j0 -1; ψ j,k , j ≥ j 0 , k = 0, 1, . . . , 2 j -1},
is an orthonormal basis of H. Thus, any function z ∈ H can then be decomposed in terms of this orthogonal basis as

z(t) = 2 j 0 -1 k=0 c j0,k φ j0,k (t) + ∞ j=j0 2 j -1 k=0 d j,k ψ j,k (t), (2)
where c j,k and d j,k are called respectively the scale and the wavelet coefficients of z at the position k of the scale j defined as

c j,k =< z, φ j,k > H , d j,k =< z, ψ j,k > H .
To efficiently calculate the WT, Mallat introduced the notion of mutiresolution analysis of H (MRA) and designed a family of fast algorithms (see Mallat (1999)).

With MRA, the first term at the right hand side of (2) can be viewed as a smooth approximation of the function z at a resolution level j 0 . The second term is the approximation error. It is composed by the aggregation of the details at scales j ≥ j 0 . These two components, approximation and details, can be viewed as a low frequency (smooth) nonstationary part and a component that keeps the time-localized details at higher scales. The distinction between the smooth part and the details is determined by the resolution j 0 , that is the scale below which the details of a signal cannot be distinguished. We will focus our attention on the finer details, i.e. on the information at the scales {j : j ≥ j 0 }.

From a practical of view, each function is usually observed on a fine time sampling grid of size N . In the sequel we will be interested in input signals of length N = 2 J for some integer J. If N is not a power of 2, one may interpolate data to the nearest J with 2 J-1 < N < 2 J . We have already seen that an advantage of the nested structure of a multiresolution analysis is that it leads to an efficient tree-structured algorithm for the decomposition of functions in V J (Mallat (1989)) for which the coefficients < z, φ J,k > H are given and that allows to derive the coefficients of the Discrete Wavelet Transform (dwt). However, when a function z is given in sampled form there is no general method for deriving the coefficients < z, φ N,k > H and one has to approximate the projection P VJ by some operator Π J in terms of the sampled values z = {z(t l ) : l = 0, . . . , N -1} of z. For regular enough wavelets, such an approximation is highly efficient (see Antoniadis (1994)) and justifies the following.

Denote by z = {z(t l ) : l = 0, . . . , N -1} the finite dimensional sample of the function z. For the particular level of granularity given by the size N of the sampling grid, one rewrites the approximation Π J of the projection P VJ of z in (2) using the truncation imposed by the 2 J points and the coarser approximation level j 0 = 0, as:

z J (t) = c 0 φ 0,0 (t) + J-1 j=0 2 j -1 k=0 d j,k ψ j,k (t), (3)

RR n°7515

where c 0 and d j,k are now denoting the empirical wavelet coefficients derived from applying the dwt on the sampled values. Hence, for a chosen, regular enough, wavelet ψ and a coarse resolution j 0 = 0, one may define the dwt operator:

W ψ : R N → R N , z → (d 0 , . . . , d J-1 , c 0 ) ,
with d j = {d j,0 , . . . , d j,2 j -1 }. Since the dwt operator is based on an L 2orthonormal basis decomposition, Parseval's theorem states that the energy of a square integrable signal is preserved under the orthogonal wavelet transform:

z 2 2 = c 2 0 + J-1 j=0 2 j -1 k=0 d 2 j,k = c 2 0 + J-1 j=0 d j 2 2 . (4)
Hence, the global energy z 2 2 of z is broken down into a few energy components. The representation ( 4) is in fact composed by the components of the discrete wavelet scalogram as defined in Arino, Morettin and Vidakovic (2004) and may be seen as the (dwt) analogue of the well-known periodogram from the spectral analysis of time series. Just as the periodogram produces an ANOVA decomposition of the energy of a signal into different Fourier frequencies, the scalogram decomposes the energy into "level components". Since N = 2 J no more than J such levels can be defined. After removing from each continuous time series slowly varying trends and eventual periodicities in time by disregarding the approximation coefficient c 0 , the scalogram components indicate at which levels of resolution the energy of the observed function is concentrated. A relatively smooth function will have most of its energy concentrated in large-scale levels, yielding a scalogram that is large for small j and small for large j. A function with a lot of high frequency oscillations will have a large portion of its energy concentrated in high resolution wavelet coefficients. Therefore the way these energies components are distributed and contribute to the global energy of a signal is the key fact that we are going to exploit to generate a handy number of features that are going to be used for clustering.

The image of the dwt operator applied on the column vector z of dimension N = 2 J may be written in matrix form as:

W = Wz,
where W is a N by N square matrix defining the dwt and satisfying W ′ W = I N (see Percival & Walden (2006, chap. 4)), and W is a column vector of length

N with W = (d 0 , d 1 , . . . , d J-1 , c 0 ) ′ ,
where W ′ denotes the transpose of W . It is easy to see that if we consider a vector x = a + bz with a, b ∈ R, then the wavelet coefficients of the dwt of x are obtained from those of the z as:

(bd 0 , bd 1 , . . . , bd J-1 , a + bc 0 ) ′ .

(5)

Absolute and relative contributions

We just have seen that dwt coefficients describe properties of functions both at various locations and at various time granularities. Each time granularity here refers to the level of detail that can be captured by dwt. This is therefore the reason of choosing the dwt as a representation scheme in our previous section to compare the shapes of curves for clustering. The energy E z = z 2 H of the time series z via decomposition (4) is equal to the sum of the energy of its wavelet coefficients distributed across scales

E z ≈ z 2 2 = c 2 0 + J-1 j=0 d j 2 2 , (6)
the approximation (denoted informally by ≈) holding because of the truncation at scale J for the wavelet expansion of z, discarding finer scales. If we consider z as the difference between two sampled curves, (6) justifies using the energy decomposition of wavelet coefficients for computing squared Euclidean distances between two series. However, when interested to see how the energy of wavelet coefficients is distributed across scales, other distance functions on dwt decompositions may be more appropriate for measuring the dissimilarity between two series.

In what follows, we define for j = 0, . . . , J -1 the absolute and relative contributions of the scale j to the global energy of the centred function respectively as

cont j = d j 2 2 , rel j = cont j J-1 j=0 cont j , ∀j = 0, . . . , J -1. (7)
We call these representations: the absolute contribution (AC) and the relative contribution (RC). We will therefore characterize each time series by the vector of its energy contributions or its relative contributions in order to define an appropriate measure of similarity that is going to be used for clustering. Note that in both of these choices of representation we leave out the eventual mean level differences of the time series since we do not make any use of the approximation term c 0 in their definition. Indeed, when within a same similarity class, the discrete time function valued processes Z are considered as stationary and therefore the scaling coefficient c 0 does not have any discriminative power, that's why we use only details after 0 in defining our distance. Now in order to compare two paths and since they are zero-mean, our distance is more relevant to measure a difference on how their energies differ across scales. In our ultimate goal to predict the future behaviour of the discrete time functional valued process Z, collecting the series of scaling coefficients for each segment Z i , i = 1, . . . , n one ends up with a real valued process of approximation coefficients. For a sequence of such univariate approximation coefficients, say {c 0,1 , . . . , c 0,n }, more or less classical time series models can be used to predict the next coefficient c 0,n+1 . Using this fact, the considered dissimilarities computed throughout the representations will be invariant under vertical shifts of the curves. Moreover, using RC implies fixing the energy of the curves to one. Hence no difference in amplitude can be detected.

3 A k-means like functional clustering procedure

We have presented a way to represent the infinite-dimensional original objects in J features that summarize the time evolution of the curves at different scales.

We will now see how we use the information that we have coded to effectively cluster it. This section starts with a brief review of some recent and highly efficient developments on feature selection and on the choice of the number of cluster.

Feature selection

Nothing warranties that all the extracted features are relevant to discover the cluster structure. Analogously with regression analysis, a feature selection step can be performed to detect the significant ones.

Feature extraction and features selection have really different aims. Whether the former creates some new information from existing objects, the latter only selects a subset of existing features. This selection reduces the computational time of the algorithm and helps to avoid an unsatisfactory and unstable clustering. Another important advantage of using a feature selection algorithm is that the reduced number of features helps to better understand the cluster output. In our case, the number of features depends on the number of sampling points of the acquired data. For N points, the number of features is J = log 2 (N ) that can be large. Moreover, since we are interested in the energy decomposition across scales, potentially several scales will not be informative for the cluster structure. Besides this, the feature selection algorithm aims to reduce (or eliminate) the presence of nonsignificant variables and a possible redundant information that could hide the cluster structure.

The absence of class labels on unsupervised learning rends particularly difficult the feature selection task. Besides, this task is intricately connected with the determination of the number of cluster. Thus, it is desirable to conduce both of them simultaneously. A recent comparative study (Steinley & Brusco (2008)) evaluate the performance of eight feature selection algorithms on simulated data sets. The compared algorithms covers model-based approaches (e.g. Law et al. (2004) which allow the simultaneous detection of groups and features) and nonparametric ones. The same authors proposed in Steinley & Brusco (2008) an algorithm that combines a variable transformation with a variable selection technique. The variable transformation introduces a variance-to-ratio weighting that looks for placing the variables on the same scale while preserving their ability to reveal cluster structures. Then, the transformed variables are used to construct an index of clusterability that serves to screen the variables that do not reveal information about the cluster structure (which is useful when working with large data sets with many masking features). Then, for the remaining variables an exhaustive evaluation of the feasible subsets of variables is done. For each subset size s, a best set of variables is obtained in terms of the largest proportion of explained variation VAF(s) from the clustering. Note that VAF decreases monotonically as a function of s. Finally, the solution of the algorithm RR n°7515

is the subset of variables that maximizes the following ratio VAF(s) -VAF(s + 1) VAF(s -1) -VAF(s)

,

where the heuristic is that for the right number of features, say s * , adding one extra variable produces a larger degradation on the VAF than the one produced when passing from s * -1 to s * .

Determination of the number of clusters

One of the most difficult task in clustering is the determination of the number of clusters K. Even if some statistical support can be given to achieve this task, usually the knowledge on the particular application helps on the choice. In the classical case, i.e. not the functional one, a lot of data-driven strategies can be defined. The first one by inspecting basically the within-cluster dissimilarity as a function of K. Many heuristics have been proposed trying to find a "kink" in the corresponding plot. A more formal argument has been proposed by Tibshirani et al. (2001) by comparing, using the gap statistic, the logarithm of the empirical within-cluster dissimilarity and the corresponding one for uniformly distributed data. Slight modifications have been proposed to the original argument, see for instance Ye (2007).

Another point of view useful to determine the number of clusters comes from model-based clustering. The idea is to fit a Gaussian mixture model to data and identify clusters as mixture components. The number of clusters is usually obtained using the well-known BIC criterion. All the above mentioned strategies seem to perform well when data do come from a mixture model but can perform poorly when the situation is more confused and fuzzy.

For determining the number of clusters James and Sugar (2003) proposed an information theoretic approach. They consider the transformed distortion curve d -p/2 K , a kind of average Mahalanobis distance between data and the set of cluster centers C as a function of K. Jumps in the associated plot allow to select sensible values for K while the largest one can be the best choice for a mixture of p-dimensional distributions with common covariance. An asymptotic analysis (as p goes to infinity) states that, when the number of clusters used is smaller than the correct number (when any), then the transformed distortion remains close to zero, before jumping suddenly and increasing linearly. Then, detecting a jump in the transformed distortion curve is equivalent to detecting the number of clusters K.

The actual procedure

We can now sum up the actual strategy for clustering that we will use in the first part of the paper in the following steps: 0. Data preprocessing. Approximate sample paths of z 1 (t), . . . , z n (t) by the truncated wavelet series (see (3)) at the scale J from sampled data z 1 , . . . , z n .

1. Feature extraction. Compute either the energetic components using absolute contribution (AC) or relative contribution (RC). If using the latter, transform the obtained vector using the logit transformation.

RR n°7515

http://www.edf.com RR n°7515


Sandra Bergonnier

Jean-Baptiste Rieunier

Stéphane Roux

Strain heterogeneities and local anisotropy in crimped glass wool

Keywords: compression test, digital image correlation, orientation analysis, strain localization

A Digital Image Correlation method is used to study a compression test carried out on crimped glass wool samples to determine the displacement field, which reveals strain heterogeneities and further localization. The structure of the wool is locally strongly anisotropic, even though the preferential orientation varies significantly in space. This structure is analyzed using an ad hoc orientation detection tool. It is shown that the strain amplitude can be correlated with texture gradients for this type of material.

Introduction

Mineral wool is a naturally textured material, well-known for its thermal and acoustic insulation properties. For some particular applications, mechanical performances are also required such as compression and tearing strengths. To evaluate the mechanical performances of the material and to quantify the effect of process parameters, mechanical tests are carried out.

To analyze tests performed on such a soft material, a contact-less method is needed. Moreover the method has to be able to deal with high strain levels. Digital Image Correlation (DIC) is particularly suited to that type of analysis. This technique is easy to use and has proven to be an efficient tool whose resolution can be extended much below the pixel accuracy (i.e., 10 -2 pixel or below in suited cases). To reach high strain levels as for mineral wool tests with a very good resolution, a multiscale approach of DIC [START_REF] Hild | Multi-Scale Displacement Field Measurements of Compressed Mineral Wool Samples by Digital Image Correlation[END_REF] can be used. It allows for large displacements between images with both a very good robustness and resolution.

Mineral wool being a highly textured material because of its natural fibrous character, it is interesting to look for a relationship between local anisotropy and mechanical strain. A tool that enables one to estimate the anisotropy of images or zones of images and the corresponding orientation field is needed. While DIC provides a strain field from mechanical tests, a texture analysis gives access to the initial orientation field of the tested sample.

The purpose of the study is to present an analysis of a compression test using both DIC and texture analysis methods. Section 2 introduces the crimping process. Section 3 deals with the DIC method. The multiscale approach and the choice of the correlation parameters are discussed. Section 4 presents the texture analysis method and its results on, first, a test case, and then on the crimped glass wool sample. Last, the mechanical test carried out on the same crimped glass wool sample is analyzed in Section 5. Compression test is performed up to -20% global strain and the relationship between the strain field and the initial texture is investigated.

Crimped glass wool

Glass wool is a cellular solid made of glass fibers whose diameter is in the range 1-5 µm, and millimetric or centimetric length. Those fibers are entangled to form a very loose mat (a few kg/m 3 ). A binder sprayed on the fibers, and cured in an oven, freezes the arrangement of the fibers and provides some elasticity to the product. As such this gives a very good insulating material, with however little mechanical performances. For applications that are more demanding on the mechanical side, the fiber mat is processed before the curing stage to produce a structure endowing the material with a stronger mechanical resistance (to compression, tearing, or shear). This step is called "crimping". It consists in breaking the naturally layered structure of the mat, through an on-line transverse and axial compression. The denser layers are buckled in the axial compression, and produce a more isotropic fiber orientation (Fig. 1), and in particular a larger population of vertical fibers, whose mechanical stiffness and strength are much larger than for the initially layered material. This underlines the key role of the product texture for its mechanical behavior, and hence a proper characterization of both the local texture, and the local strain (and their interdependence), is required to optimize the product. This motivates the study detailed in Section 5.

Displacement measurement by digital image correlation

DIC is based on the use of an image matching algorithm (i.e., correlation product), which can be performed either in the physical space (Sutton et al., 1983;[START_REF] Chu | Applications of Digital-Image-Correlation Techniques to Experimental Mechanics[END_REF] or in Fourier space [START_REF] Chen | Digital Speckle-Displacement Measurement Using a Complex Spectrum Method[END_REF][START_REF] Berthaud | Méthodes optiques et acoustiques de mesures des caractéristiques mécaniques[END_REF]. In the present case, Fourier computations are used. The in-plane displacement field is evaluated with the cross-correlation of an interrogation window in the picture of the deformed surface with respect to a reference one. To determine the displacement field of one image of the deformed surface with respect to a reference image, one considers a sub-image that will be referred to as a Zone Of Interest (ZOI). The aim of the correlation method is to match the ZOIs in the two images (Fig. 2). The displacement of a ZOI with respect to the other one is a two-dimensional shift of an intensity signal digitized by a CCD camera. In practice, two images are considered. The first one, referred to as 'reference image' and the second one, called 'deformed image'. The reference image, schematically shown in Fig. 3, is meshed by ZOIs that form the region of interest (ROI). The sub-ROI of size 2 P × 2 Q , is extracted from the ROI in the reference image with the largest values of P and Q. First a standard approach can be used [START_REF] Périé | Analysis of a Multiaxial Test on a C/C Composite by Using Digital Image Correlation and a Damage Model[END_REF], which is wellsuited for small strain perturbations. But if displacements are too large, the information in the ZOI of the reference image will be partially lost in the ZOI of the deformed image, i.e., there is not enough information to determine a reliable correction. To use the method in the domain of large strains, as for the case of mineral wool, a large ZOI size should be used, whose main drawback is to average displacements over a larger surface. Yet the aim of a correla- tion calculation is to determine displacement fields at the scale of texture heterogeneities. The maximum measurable displacement for the algorithm to converge to the proper solution has then to be increased. This is achieved by using a multiscale algorithm.

MULTISCALE ALGORITHM

The aim of the multiscale algorithm is to improve the spatial resolution (i.e., to decrease the ZOI size), yet allowing for a large displacement between two consecutive images. The previous approach is recursively applied on images at different scales. First the average strains are deduced from an analysis with images where finer details have been erased (i.e., super-pictures, as shown in Fig. 4). To make the evaluation more robust, the interpolated displacements are applied to the analyzed ZOI and re-evaluated until the maximum difference between two iterations is smaller than a chosen value. Details can then be restored in the image (i.e., resolving finer scales) and the displacement interpolation can be used as a first evaluation of the field at the new scale. This procedure is iterated until the scale corresponds to the initial ROI (with all its details). During the last step, the displacement is no longer interpolated and fluctuations, whose amplitudes are greater than the ones that can be reached with the standard approach, can be analyzed.

Figure 4 shows the construction of images at different scales. Scale no. 0 corresponds to the considered image and scale no. 1 to the region of interest. Scale no. 2, extracted from the ROI, is determined by the largest values for P and Q so that the corresponding size of the sub-ROI is 2 P × 2 Q (Fig. 3). From scale no. 3 on, each transition is characterized by the definition of super-pixels. The latter are defined recursively from one scale to the next by averaging the gray levels of 2 × 2 neighboring pixels. This procedure is carried out until the minimum size of the sub-image is equal to 128 pixels. The analysis starts at, say, scale no. 4 and ends at scale no. 1. This allows one to obtain a maximum measurable strain four times larger than the standard approach even if the maximum measurable displacement is still limited [START_REF] Hild | Multi-Scale Displacement Field Measurements of Compressed Mineral Wool Samples by Digital Image Correlation[END_REF], as in the standard approach, only to the coarser scale (here scale no. 4). On artificially deformed images, strains of the order of ±30% can be measured [START_REF] Hild | Multi-Scale Displacement Field Measurements of Compressed Mineral Wool Samples by Digital Image Correlation[END_REF][START_REF] Chevalier | Tools for Multiaxial Validation of Behavior Laws Chosen for Modeling Hyper-Elasticity of Rubber-Like Materials[END_REF] between two pictures. When a sequence of more than 2 images is analyzed and when the reference image is updated (i.e., it corresponds to the deformed image of the previous step), there is no real limitation, apart from the fact that errors are cumulated.

CHOICE OF THE CORRELATION PARAMETERS

The aim of this part is to estimate the uncertainty of the correlation algorithm used herein. The displacement and strain uncertainties are estimated when the correlation parameters are modified. To have independent estimates of the displacements in the following analysis, the shift δ is equal to the ZOI size (i.e., δ = 2 n pixels).

First, the displacement uncertainty is estimated. In the present case, a constant displacement of 0.5 pixel is applied artificially by using the shift/modulation property of Fourier transforms [START_REF] Périé | Analysis of a Multiaxial Test on a C/C Composite by Using Digital Image Correlation and a Damage Model[END_REF]. The average difference gives an evaluation of the error, and the standard deviation that of the corresponding uncertainty. These properties are associated to the analyzed texture (Fig. 1). It can be shown that the displacement uncertainty is half that reached when the prescribed displacement is equal to 0.5 pixel [START_REF] Bergonnier | Digital image correlation used for mechanical tests on crimped glass wool samples[END_REF]. Figure 5 shows the displacement uncertainty for different ZOI sizes. The larger the ZOI, the smaller the error. However, if the ZOI size increases, the spatial uncertainty is degraded, i.e., less independent measurement points are obtained.

Second, the strain uncertainty is determined by using the same picture. From the displacement field, in-plane strains are evaluated by using a classical finite differences scheme so that the spatial uncertainty (i.e., the gauge length) is equal to 2 × δ, i.e., 2 n+1 pixels. In all cases, the mean strain is less than 7 × 10 -5 . The strain uncertainty is more sensitive to δ (here equal to the ZOI size) than the displacement uncertainty to the ZOI size (Fig. 5).

A good compromise for the present study is given by a ZOI size of 16 pixels (i.e., n = 4) for which a displacement uncertainty of the order of 0.03 pixel is achieved. In terms of strains, the corresponding uncertainty is 3×10 -3 with a spatial uncertainty of 32 pixels. It should be remembered that all the results discussed so far were obtained with the same reference picture. When large strains or severe texture variations occur [START_REF] Hild | Multi-Scale Displacement Field Measurements of Compressed Mineral Wool Samples by Digital Image Correlation[END_REF], an updating procedure is preferred. Consequently the uncertainty gradually degrades as the number of analyzed pictures increases.

Local anisotropy analysis

Detecting the anisotropy of a domain, contained in images, involves an evaluation of the change of gray levels in different directions, the direction where the variation is the least being the principal axis of anisotropy. Since changes in gray levels are to be considered, most methods [START_REF] Rao | A Taxonomy for Texture Description and Identification[END_REF][START_REF] Germain | Multiscale estimation of vector field application to texture characterization[END_REF] are based on the calculation of gradient fields and implicitly assume a smoothness of the image. The method that is presented here does not require any smoothness assumption of the image. This may be an important improvement for images where the anisotropy is provided by very high frequency modes (i.e., fine texture).

Let f (x, y) represent a gray-scale image. The variation of gray levels suggests the use of an L 2 -norm of the finite difference along the e direction

D e f =   (x,y) (f (x + e x , y + e y ) -f (x, y)) 2   1/2 (1)
where e x , e y are the components of the vector e. This L 2 -norm can easily be computed for any direction in Fourier space. The total power is preserved in the Fourier transform, and thus

D e f 2 = k (k.e) 2 | f (k)| 2 (2)
Therefore, introducing the anisotropy tensor T defined as

T = k (k ⊗ k)| f (k)| 2 (3) one observes that D e f 2 = T : (e ⊗ e) (4)
where ':' is the tensorial contraction with respect to two indices. The anisotropy axis appears as the minor eigendirection of the anisotropy tensor T. The anisotropy amplitude can be estimated conveniently through dev(T)/tr(T), where 'tr' designates the and 'dev' the norm of the deviator T-tr(T)I/2.

APPLICATION TO A TEST IMAGE

The method presented above is tested on concentric circles that represent an absolute reference because the true orientation is known at each position in the image. The image size is 512 × 512 pixels, and the size of the zones of interest 16 × 16 pixels. These zones are centered on a regular square grid and the distance between the center of two consecutive zones is 16 pixels as for the DIC analysis (Fig. 3). In Fig. 6(a), where only one quarter of the image is shown, a first visual inspection gives a rough evaluation of the adequation between real orientations and estimates. Figure 6(b) shows the differences between the determined and true orientations as a function of the angle. The 0 • angle corresponds to a horizontal orientation. The RMS error or resolution is 0.63 • .

APPLICATION TO CRIMPED GLASS WOOL

The is then applied to the determination of the texture orientation of a crimped glass wool sample. A special strategy has been used to limit the influence of edge effects (Bergonnier et al., 2005). For such an image (Fig. 1), is to be noticed that there is a very fine texture, with a very difficult partition between noise and information carrying signal. Moreover, if anisotropy is clearly marked, the coherence length is rather short, and the pattern appears as circumvoluted. This is thus a very severe test of the proposed method. In Fig. 7, the texture orientation corresponding to the crimped glass wool sample (Fig. 1) is presented. A discontinuity of the orientation between the upper and the lower parts of the sample can be observed, with a 50 • orientation for the upper part and a -50 • orientation for the lower part. The fact that the texture is very fine can also be observed in Fig. 7 where the variations in the orientation field can be important from one ZOI to its neighbors. If the local texture of the sample has to be linked with the mechanical response of the sample to a compression test, the DIC method will enable one to compare the strain field and the initial orientation field by using the same calculation points.

Experiments: Uniaxial compression test on crimped glass wool

The compression test reported herein is displacement-controlled (or equivalently, the global strain ǫ is controlled; it is defined as the ratio of the displacement between the two platens divided by the initial height of the sample). A uniaxial compression (in the 1-direction) up to -20% global strain is applied on parallelepipedic samples, whose dimensions are 200 × 200 × 80 mm 3 , 80 mm being the height of the sample (Fig. 1). For the test, images of the sample illuminated with direct light have been taken for different global strain levels (with -1% increments in compression).

result of the test is shown in Fig. 8. The "plastic" regime suggests the existence of heterogeneous strains. This implies to study the strain field of the material with a good spatial resolution. In Fig. 9, initial texture orientation and nominal strain fields are shown. For the latter, the analyzed kinematic quantity is evaluated for a reference given prior to the beginning of the experiment. The ZOIs used for the calculation correspond to the ZOIs used for the local anisotropy analysis. The calculation, performed with ZOIs of size 16 × 16 pixels with a 16-pixel shift δ, corresponds to the component ǫ 11 for different global strain levels with an adaptive scale, namely, -5% for Fig. 9(b), -7% for Fig. 9(c), -10% for Fig. 9(d) and -15% for Fig. 9(e), respectively. In Fig. 8, the corresponding points are located in the elastic domain, at the end of the elastic domain, at the beginning of the "plastic" domain and where the "plastic" domain is well-defined. For the calculation of the strain fields the same filtering as for the anisotropy analysis has been applied. Within the elastic domain (Fig. 9(b-c)) a first strain localization zone can be identified that seems to traverse the sample. For a -10% global strain level (Fig. 9(d)), strain heterogeneities can well be observed at the surface of the material and a well-marked localized zone appears that correspond to the first zone identified in Fig. 9(b-c). For -15% global strain level, strain fluctuations appear around the previously observed localized zone.

The strain localization zone corresponds to a band that traverses the surface of the sample. If the orientation map (Fig. 9(a)) is compared to the strain fields, one can observe that the discontinuity in the orientation field is located where strain localization band occurs. The well-marked discontinuity zone on the left part of the material corresponds to where the first strain localization zone can be observed at early stages within the elastic domain (Fig. 9(b-c)).

Conclusion

First, a multiscale approach of the displacement and strain fields identification has been presented. The algorithm increases both the robustness and the maximum measurable displacement and strain between two consecutive images. This technique is well-suited to analyzing compression tests carried out on crimped glass wool samples.

Second, an image analysis tool has been introduced to determine the local anisotropy of an image. This approach gives access to the orientation field of an image. Validated on an artificial case, it was then applied on a crimped glass wool sample picture.

The results of a compression test on a crimped glass wool sample are studied. The DIC technique enables one to identify strain heterogeneity zones, and then strain localization bands. If both methods, DIC and local anisotropy analysis are coupled, the discontinuity zone in orientation can be linked to the localization bands in the strain field.

The correlation technique applied to crimped glass wool shows the change of macrostructure during the test. A further development is to characterize the change in the local orientation by the presented method and to relate it to strain localization.

Figure 1 .

 1 Figure 1. Lateral view of a crimped glass wool product (1 pixel ⇔ 0.15 mm).

Figure 2 .

 2 Figure 2. Zone of Interest (ZOI) in reference (left) and "deformed" (right) images. Associated displacement vector.

Figure 3 .

 3 Figure 3. Schematic of the reference image with the correlation parameters (ZOI area: 2 n × 2 n pixels and shift δ).

Figure 4 .

 4 Figure 4. Different scales used in the multi-scale approach for a glass wool sample.

Figure 5 .

 5 Figure 5. Displacement and strain uncertainties for different correlation parameters (ZOI size = 2 n pixels, shift δ = 2 n pixels).

Figure 6 .

 6 Figure 6. (a) Results of the anisotropy analysis on a test image, (b) Differences between calculated and true orientations.

Figure 7 .

 7 Figure 7. Texture orientation (in degrees) of the crimped glass wool sample before the mechanical test.

Figure 8 .

 8 Fig. 9(b)

Figure 9 .

 9 Figure 9. Initial texture orientation and cumulative ǫ11-field for different levels of the global strain. Correlation parameters: size of ZOI = 16 × 16 pixels, δ = 16 pixels. Scale: pixel ⇔ 0.15 mm.

  Sutton, M.A., Wolters, W.J., Peters, W.H., Ranson, W.F. and McNeill, S.R.: Determination of Displacements Using an Improved Digital Correlation Method. Im. Vis. Comp., 1(3): 133-139, 1983.

Bergonnier.tex; 10/01/2005; 17:49; p.12

Bergonnier.tex; 10/01/2005; 17:49; p.2

Bergonnier.tex; 10/01/2005; 17:49; p.3

Bergonnier.tex; 10/01/2005; 17:49; p.4

Bergonnier.tex; 10/01/2005; 17:49; p.5

Bergonnier.tex; 10/01/2005; 17:49; p.6

Bergonnier.tex; 10/01/2005; 17:49; p.7

Bergonnier.tex; 10/01/2005; 17:49; p.8 Strain heterogeneities and local anisotropy in crimped glass wool

Bergonnier.tex; 10/01/2005; 17:49; p.9

Bergonnier.tex; 10/01/2005; 17:49; p.10

Bergonnier.tex; 10/01/2005; 17:49; p.11

Acknowledgements

This work has been supported by Isover Saint-Gobain and by CNRS in a project entitled "Analyses multi-échelles de champs de déformation par traitement d'image : vers l'identification de champs de propriétés mécaniques." The authors acknowledge useful discussions with M. Baudequin, O. Gaume and Machelart.


Evangelos Maria Dimopoulou

Prof Evangelos Terpos

Ersi Voskaridou
email: ersi_voskaridou@yahoo.com

Eleni Plata

Marousa Douskou

Anastasia Sioni

Efrosini Mpoutou

Dimitrios Christoulas

MD Evangelos Terpos

Deferasirox effectively decreases iron burden in patients with double heterozygous HbS/-thalassemia

Keywords: double heterozygous HbS/beta-thalassemia, sickle cell syndromes, deferasirox, iron overload, liver, heart double heterozygous HbS/-thalassemia, sickle cell syndromes, deferasirox, iron overload, liver, heart

Deferasirox effectively decreases iron burden in patients with double heterozygous HbS/β-thalassemia

Ersi Voskaridou, Eleni Plata, Marousa Douskou, Anastasia Sioni, Efrosini Mpoutou, Dimitrios Christoulas,

are using these values and related cut-offs for the assessment of liver siderosis." In answering this comment the authors changed the text and included a reference explaining the use of R2 and R2* values. Actually the cited paper by J. Wood et al. accurately gives information on the hepatic iron concentration (HIC) in [mg/g d.w.] in relation to the corresponding R2 or R2* values. Such information on HIC corresponding to liver T2 and T2* values which were used in the present study would be desirable. What HIC corresponds e.g. to the inclusion criterion of T2* of 19ms? Answer: We have included the respective HIC in the inclusion criteria (page 4, lines 5-6) as well as in the Table 1 (4th line), according to reviewer suggestions. Comment 2. (see previous comment 6) : In "Methods" it is stated that: "The starting deferasirox dose was 10 or 20 mg/kg/day depending on baseline iron burden..." In "Results", it is written that "patients were treated with deferasirox at 10 and 20mg/kg/day? based on the number of transfusions?")? How was the "baseline iron burden defined" (according to "Results" by the number of transfusions) and what was the cut-off number of blood transfusions to decide for 20mg/kg/d? Answer: We thank the reviewer for this comment and sorry that we had not realized it in our previous revision. We have clarified now in the "study design" section that "The starting deferasirox dose was 10 or 20 mg/kg/day depending on previous transfusion requirements: if the patient had received 20-30 red blood cell (RBC) units the dose of 10 mg/kg/day was used; if the patients had received more than 30 RBC units the dose of 20 mg/kg/day was used." [page 4, lines 15-17].

Dear Professor Kulozik

We would like to thank you and the reviewers for the re-evaluation of our paper entitled "Deferasirox effectively decreases iron burden in patients with double heterozygous HbS/-thalassemia". Please find enclosed the second revised edition of this paper based on the reviewer's comments. All changes have been bolded into the text.

Reviewer #1:

Comment 1. The manuscript was thoroughly revised and has remarkably improved. However, concerning two of the previous comments were incompletely answered. A more detailed answer will increase the comprehensibility of the corresponding parts of the paper. Certainly, it will be very easy for authors to add the missing information: As stated in the previous comment 5, for liver MRI, it would be helpful for the reader to know how the T2star values correspond to liver iron concentration given as mg/g or micromole/g dry or wet weight since at present the majority of institutions and colleagues are using these values and related cut-offs for the assessment of : In "Methods" it is stated that: "The starting deferasirox dose was 10 or 20 mg/kg/day depending on baseline iron burden..." In "Results", it is written that "patients were treated with deferasirox at 10 and 20mg/kg/day? based on the number of transfusions?")? How was the "baseline iron burden defined" (according to "Results" by the number of transfusions) and what was the cut-off number of blood transfusions to decide for 20mg/kg/d? Answer: We thank the reviewer for this comment and sorry that we had not realized it in our previous revision. We have clarified now in the "study design" section that "The starting deferasirox dose was 10 or 20 mg/kg/day depending on previous transfusion requirements: if the patient had received 20-30 red blood cell (RBC) units the dose of 10 mg/kg/day was used; if the patients had received more than 30 RBC units the dose of 20 mg/kg/day was used." [page 4, lines 15-17].

We think that we have addressed the reviewer comments adequately and that now our paper will be accepted for publication in Annals of Hematology.

We look forward to receiving your final decision on our paper.

Chairman of Thalassemia Unit

Introduction

The clinical severity of sickle-cell disease (SCD) varies from a very mild disorder to severe organ failure with increased mortality. Similarly, the double heterozygosity of SCD and betathalassemia (HbS/thal) produces a variety of clinical pictures with those of SCD to be dominant [START_REF] Charache | Effect of hydroxyurea on the frequency of painful crises in sickle cell anemia. Investigators of the Multicenter Study of Hydroxyurea in Sickle Cell Anemia[END_REF][START_REF] Steinberg | Effect of hydroxyurea on mortality and morbidity in adult sickle cell anemia: risks and benefits up to 9 years of treatment[END_REF][START_REF] Voskaridou | The effect of prolonged administration of hydroxyurea on morbidity and mortality in adult patients with sickle-cell syndromes: results of a 17-year, single center trial (LaSHS)[END_REF]. Despite adequate therapy, the majority of patients with HbS/thal have received repeated blood transfusions by adulthood [START_REF] Wanko | Transfusion management in sickle cell disease[END_REF][START_REF] Adams | Prevention of a first stroke by transfusions in children with sickle cell anemia and abnormal results on transcranial Doppler ultrasonography[END_REF]. An imbalance between chronic inflammation which increases hepcidin expression, chronic hemolysis, accumulation of transfusional iron and increased iron absorption due to ineffective erythropoiesis leads to disruption of iron metabolism and possibly to iron overload in several HbS/thal patients. The increase in longevity during the recent years [START_REF] Steinberg | Effect of hydroxyurea on mortality and morbidity in adult sickle cell anemia: risks and benefits up to 9 years of treatment[END_REF][START_REF] Platt | Mortality in sickle cell disease (1994) Life expectancy and risk factors for early death[END_REF] has revealed clinical evidence of iron overload in several patients with SCD and HbS/thal, especially in the liver; thus iron chelation has a definite indication in these cases. Deferasirox (Exjade) is a novel iron chelator that has been licensed for the management of transfusional iron overload. There are several reports for the efficacy of deferasirox in thalassemia, SCD and other transfusion-dependent anemias [START_REF] Angelucci | Italian Society of Hematology practice guidelines for the management of iron overload in thalassemia major and related disorders[END_REF][START_REF] Voskaridou | Treatment with deferasirox (Exjade®) effectively decreases iron burden in patients with thalassaemia intermedia: results of a pilot study[END_REF][START_REF] Cappellini | Tailoring iron chelation by iron intake and serum ferritin: prospective EPIC study of deferasirox in 1744 patients with transfusion-dependent anemias[END_REF][START_REF] Vichinsky | Deferasirox in Sickle Cell Investigators (2007) A randomised comparison of deferasirox versus deferoxamine for the treatment of transfusional iron overload in sickle cell disease[END_REF] but there are very limited data for its role in HbS/thal. The aim of this study was to investigate the efficacy and safety of deferasirox on iron-overloaded patients with HbS/thal.

Patients and methods

Inclusion and exclusion criteria: Adult patients (≥18 years of age) with HbS/thal , who had received sporadic red blood cell transfusions (≤20 units of RBCs in their lifetime) and had confirmed liver or cardiac iron overload (serum ferritin levels of ≥1000 ng/mL and one of the following: liver MRI T2* of <19 ms, which corresponds to a liver iron concentration (LIC) of 1.5 mg/g dry weight or cardiac T2* of <28 ms) were eligible for entrance into this study.

Patients were excluded if they were pregnant, showed signs of hepatic failure (transaminase levels of >500 U/L) or renal failure (creatinine clearance <60 mL/min) or had a left ventricular ejection fraction (LVEF) of <50%. Regarding the cut-off value of MRI T2* of the heart for inclusion in the study, we used the level of 28 ms instead of 20 ms which is used in the majority of the studies, as in another study of our group none of the patients with SCD had a heart T2* MRI value of <28 ms [START_REF] Voskaridou | Angiogenic molecules and inflammatory cytokines in patients with thalassemia major and double heterozygous HbS/beta-thalassemia; the impact of deferasirox[END_REF].

Study design:

This was an open-label, prospective, single-centre, phase II study on the role of deferasirox in iron overload of HbS/thal, conducted over 1-year. The starting deferasirox dose was 10 or 20 mg/kg/day depending on previous transfusion requirements: if the patient had received 20-30 red blood cell (RBC) units the dose of 10 mg/kg/day was used; if the patients had received more than 30 RBC units the dose of 20 mg/kg/day was used. Dose adjustments were permitted after 3 months (in increments of 5 mg/kg/day every 3 months as required) based on serum ferritin trends (levels of ≥1000 ng/mL on 2 sequential visits or levels of >2500 ng/mL without decreasing trend). Deferasirox dosage was reduced if there were elevated levels of creatinine, urinary protein:creatinine ratio and transaminases, and in response to adverse events. Patients provided written, informed consent before entering the study. The study was conducted in accordance with Good Clinical Practice guidelines and the Declaration of Helsinki. Hemoglobin (Hb) levels, serum creatinine, alanine (ALT) and aspartate aminotransferase (AST) and 24-hour proteinurea (Prot 24h) were measured before and every month during deferasirox treatment, using standard methodology. Serum cystatin C was also measured at the same timepoints on a Behring Nephelometer-II analyser using a latex particle-enhanced nephelometric immunoassay (Dade Behring, Liederbach, Germany), as previously described [13].

Efficacy and safety assessments:

MRI evaluation: All patients were examined using the same 1.5 Tesla apparatus (Philips 1.5 T Gyroscan Intera; the Netherlands). Morphological evaluation and gross pathological assessment of liver was based on high resolution, T2-weighted, turbo spin echo images (TR/TE/matrix: 1600ms/100ms/400). For spin-spin relaxation time (T2) measurements, images of the liver were obtained using a 2D, 16-echo, Carr-Purcell-Meiboom-Gill spin echo sequence (TR/TE/matrix : 1000ms/5ms/256) T2-maps, so that images in which pixel signal intensity is related to the T2 of the corresponding voxel, can be automatically derived by the imager's software. These maps were used for hepatic parenchyma T2 determination within selected regions of interest. For the T2 measurements of the heart, a short axis midventricular slice was acquired. Respiratory and cardiac triggering were used for the heart measurements. T2* MRI of the heart and the liver were performed as previously described [START_REF] Voskaridou | Angiogenic molecules and inflammatory cytokines in patients with thalassemia major and double heterozygous HbS/beta-thalassemia; the impact of deferasirox[END_REF]. LVEF measurements were also performed using a phased array cardiac coil and a multislice, multi-phase ECG, gated bFFE sequence. Normal ranges for the T values of liver and heart were established in 10 age-matched volunteers (4 males and 6 females).

Statistical analysis:

The Wilcoxon signed rank test was used to evaluate differences between baseline and values of the studied parameters at the various time points. The correlation between changes of various biochemical parameters and BMD was evaluated with the Spearman's (r s ) correlation coefficient. All p values are two sided, the level of significance is <0.05 and confidence intervals refer to 95% boundaries.

Results

Thirty-one adult patients with HbS/thal (14M/17F) who are regularly followed-up in the Thalassemia Center of Laikon Hospital (Athens, Greece) were entered into this study. The median age of all patients was 41 years (range: 27-55 years). Twenty-four patients had HbS/  thal, while 7 patients had HbS/  thal. Twenty-three and 8 patients were initially treated with deferasirox at 10 and 20 mg/kg/day, respectively. Deferasirox produced a significant reduction of mean serum ferritin levels after 6 and 12 months of treatment, while the mean baseline liver T2 and T2* significantly increased following 12 months of therapy (Table 1 & Figure 1). Mean cardiac T2* and LVEF were normal at baseline and did not change after 12 months of treatment (normal values of T2* were >28 ms and of LVEF >65%). A strong correlation was observed at baseline between liver MRI T2* and serum ferritin (r=0.61, P<0.01), but there was only weak correlation between percentage changes of serum ferritin and liver MRI T2* after 12 months of therapy (r=0.233, P=0.09).

The reduction of ferritin levels after 12 months of deferasirox administration (approximately 50% of the mean values) was accompanied by a comparable reduction in the mean levels of ALT and AST (40% and 30%, respectively; Table 1). There were no significant correlations between serum ferritin and heart MRI T2* or between MRI T2* of heart and liver at baseline or after the 12-month duration of the study.

Regarding renal function, there were no significant changes in mean serum creatinine; all observed increases were within the normal range of our laboratory. On the contrary, cystatin-C significantly increased after 6 months of treatment and remained increased post-12 month of deferasirox. Hb or Prot 24h levels did not change after 6 and 12 months of treatment.

The side-effects of deferasirox were minimal and easily manageable. Eight patients (25%) referred gastrointestinal disturbances: 5 had mild nausea and 4 had diarrhea (one had both).

All were reported within the first month, they were of grade I or II, and they did not reoccur. In all patients the compliance was excellent.

There were no dose adjustments of deferasirox in this study.

Discussion

Iron overload may play an important role in increasing the morbidity and mortality of patients with HbS/thal and iron chelation has a definite indication in several cases. Chelation guidelines for sickle-cell syndromes are similar to those for other iron overload situations as thalassemia major and multi-transfused patients [14,15]. Until recently, desferrioxamine and deferiprone were the only active iron-chelating agent approved for clinical use, mainly in the thalassemia setting [16,17]. Recent advances have led to the development of new oral iron chelators, such as deferasirox, which has changed the clinical practice in several transfusiondependent anemias with iron overload. Vichinsky et al suggest that treatment adherence with deferasirox may be better than that of desferrioxamine and thus it should lead to improved long-term outcomes [18]. Deferasirox mobilizes iron stores by binding selectively to the ferric form of iron. It has been shown to be as effective as desferrioxamine with favorable safety profile in patients with different hemoglobinopathies [START_REF] Angelucci | Italian Society of Hematology practice guidelines for the management of iron overload in thalassemia major and related disorders[END_REF][START_REF] Voskaridou | Treatment with deferasirox (Exjade®) effectively decreases iron burden in patients with thalassaemia intermedia: results of a pilot study[END_REF][START_REF] Cappellini | Tailoring iron chelation by iron intake and serum ferritin: prospective EPIC study of deferasirox in 1744 patients with transfusion-dependent anemias[END_REF][START_REF] Vichinsky | Deferasirox in Sickle Cell Investigators (2007) A randomised comparison of deferasirox versus deferoxamine for the treatment of transfusional iron overload in sickle cell disease[END_REF]. However, there is very limited information in the literature for the effect of deferasirox on HbS/thal patients.

In our study deferasirox was proven very effective on reducing liver iron and serum ferritin in the vast majority of HbS/thal patients. Serum ferritin is a marker of iron overload but also an acute phase protein. Its reduction post-deferasirox may reflect a beneficial effect of deferasirox on reducing iron overload, a possible anti-inflammatory effect of deferasirox or both phenomena. Our group has recently presented that deferasirox reduces the levels of tumor necrosis-alpha (a major inflammatory cytokine) in patients with thalassemia major and HbS/-thal [19].

Iron overload contributes to the severity of SCD aggravating the already compromised liver and increasing the probability of cardiac complications. Our group has previously shown that deferiprone is able to reduce liver iron overload in sickle-cell syndromes [START_REF] Voskaridou | Deferiprone as an oral iron chelator in sickle cell disease[END_REF]. In the present study, which according to our knowledge is the first where deferasirox has been given in patients with HbS/thal, we also showed that liver T2 and T2* significantly increased following 12 months of deferasirox therapy. Additionally, the reductions in liver enzymes (ALT and AST) post-deferasirox administration further support a possible beneficial effect of deferasirox on liver function. Liver disease does occur in patients with sickle-cell syndromes not only as a result of chronic obstruction by sickled red cells but also because of fibrosis and cirrhosis induced by iron deposition [START_REF] Banarjee | Sickle cell Hepatopathy[END_REF]. Thus, in these patients, liver is considered as the main target organ and deferasirox seems to improve its function by removing the excess iron.

On the contrary, heart MRI T2* failed to display any differences between baseline and 12month values. However, we have to mention that all baseline and post-therapy values were within normal values. Furthermore, heart MRI T2* showed no correlations with either the serum ferritin levels nor with liver MRI T2*. These results are in agreement with those of Inati et al who reported that no one of 23 evaluated SCD patients had evidence of cardiac iron overload although they had significant transfusion burden, systemic and hepatic iron overload. Moreover, the authors found no significant correlations between cardiac MRI T2* values and biochemical variables [START_REF] Inati | Absence of cardiac siderosis by MRI T2* despite transfusion burden, hepatic and serum iron overload in Lebanese patients with sickle cell disease[END_REF].

Preclinical studies of deferasirox indicated that the kidney is a potential target organ of toxicity. Patients with SCD have often normal serum creatinine values, although their renal function is abnormal. In this study, we evaluated a reliable and sensitive marker of glomerular filtration rate, named cystatine-C [START_REF] Coll | Serum cystatin C as a new marker for noninvasive estimation of glomerular filtration rate and as a marker for early renal impairment[END_REF], which has been found to be increased in patients with sickle-cell syndromes, including HbS/thal [START_REF] Voskaridou | Early markers of renal dysfunction in patients with sickle cell/beta-thalassemia[END_REF]. Deferasirox increased serum cystatin-C at 6 and 12 months of therapy, but there was no further increase between the 6 th and 12 th month of treatment. Such an increase was not observed in patients with thalassemia intermedia who were treated with deferasirox in a similar fashion in our center [START_REF] Voskaridou | Treatment with deferasirox (Exjade®) effectively decreases iron burden in patients with thalassaemia intermedia: results of a pilot study[END_REF]. Vichinsky et al reported in a randomized study of 195 adult and pediatric patients with SCD that deferasirox produced mild non-progressive increases in serum creatinine of some patients [START_REF] Vichinsky | Deferasirox in Sickle Cell Investigators (2007) A randomised comparison of deferasirox versus deferoxamine for the treatment of transfusional iron overload in sickle cell disease[END_REF]. In our study, we did not notice any increases in serum creatinine values or any changes in 24-hour proteinuria. However, our results on cystatin-C levels suggest that a very close monitor of renal function is needed when we use deferasirox in HbS/thal. We conclude that our data indicate that over 12 months, deferasirox significantly reduced liver iron burden and serum ferritin levels in iron-overloaded patients with HbS/thal. Our study suggests that deferasirox provides effective iron chelation therapy of the liver in this cohort of HbS/thal patients with minimal toxicity.

  liver siderosis." In answering this comment the authors changed the text and included a reference explaining the use of R2 and R2* values. Actually the cited paper by J. Wood et al. accurately gives information on the hepatic iron concentration (HIC) in [mg/g d.w.] in relation to the corresponding R2 or R2* values. Such information on HIC corresponding to liver T2 and T2* values which were used in the present study would be desirable. What HIC corresponds e.g. to the inclusion criterion of T2* of 19ms? Answer: We have included the respective HIC in the inclusion criteria (page 4, lines 5-6) as well as in the

  The primary efficacy endpoint was the improvement in magnetic resonance imaging (MRI) T2* values of the liver and the heart from baseline to 12month values post-deferasirox therapy. The secondary efficacy endpoints included evaluation of changes in serum ferritin and liver function, and assessment of incidence, type and severity of adverse events relating to biochemical changes and patient disposition. We have used liver MRI T2* change as a primary efficacy end-point because it reflects accurately the LIC. Indeed a combination of the relaxation rates R2 (1/T2) and R2* (1/T2*) measured by MRI can accurately measure hepatic iron concentration throughout the clinically relevant range of LIC with appropriate MRI acquisition techniques [12].

  overload in patients with beta thalassaemia and sickle cell disease. Br J Haematol 126:736-742 12. Wood JC, Enriquez C, Ghugre N, Tyzka JM, Carson S, Nelson MD, Coates TD (2005) MRI R2 and R2* mapping accurately estimates hepatic iron concentration in transfusiondependent thalassemia and sickle cell disease patients. Blood 106:1460-1465 13. Terpos E, Katodritou E, Tsiftsakis E, Kastritis E, Christoulas D, Pouli A, Michalis E, Verrou E, Anargyrou K, Tsionos K, Dimopoulos MA, Zervas K; Greek Myeloma Study Group (2009) Cystatin-C is an independent prognostic factor for survival in multiple myeloma and is reduced by bortezomib administration. Haematologica 94:372-379 14. Porter JB, Huehns ER (1987) Transfusion and exchange transfusions in sickle cell anemias, with particular reference to iron metabolism. Acta Haematol 78:198-205 15. Charache S (1981) Treatment of sickle cell anemia. Annu Rev Med 32:195-206 16. Pakbaz Z, Fischer R, Treadwell M, Yamashita R, Fung EB, Calvelli L, Quirolo K, Foote D, Harmatz P, Vichinsky EP (2005) A simple model to assess and improve adherence to iron chelation therapy with deferoxamine in patients with thalassemia. Ann N Y Acad Sci 1054:486-491 17. Cohen AR, Galanello R, Piga A, De Sanctis V, Tricta F (2003) Safety and effectiveness of long-term therapy with the oral iron chelator deferiprone. Blood 102:1583-1587 18. Vichinsky E, Pakbaz Z, Onyekwere O, Porter J, Swerdlow P, Coates T, Lane P, Files B, Mueller BU, Coïc L, Forni GL, Fischer R, Marks P, Rofail D, Abetz L, Baladi JF (2008) Patient-reported outcomes of deferasirox (Exjade, ICL670) versus deferoxamine in sickle cell disease patients with transfusional hemosiderosis. Substudy of a randomized openlabel phase II trial. Acta Haematol 119:133-141

Figure 1 .

 1 Figure 1.

Table 1 (

 1

	4 th line), according to reviewer
	suggestions.
	Comment 2. (see previous comment 6)

Table 1 .

 1 All evaluated parameters at baseline, after 6 and 12 months of deferasirox treatment (mean±SD)

		Baseline	6 months	12 months
	Serum ferritin, ng/mL	1989 ± 923	1620 ± 894	1008 ± 776
			(P=0.05)	(P<0.001)
	Liver T2, ms	20.4 ± 6.0	NM	28.2 ± 8.3
				(P=0.001)
	Liver T2*, ms	3.9 ± 3.2	NM	5.8 ± 3.1
				(P=0.01)
	HIC, mg/g dry weight	6.9 ± 8.7	NM	4.7 ± 8.8
				(P=0.01)
	Cardiac T2*, ms	35.4 ± 8.3	NM	37.8 ± 6.4
				(P=0.481)
	LVEF, %	64.5 ± 7.4	NM	66.0 ± 6.8
				(P=0.376)
	AST, U/L	57.7 ± 26.8	47.9 ± 29.4	40.5 ± 22.6
			(P=0.06)	(P<0.01)
	ALT, U/L	50.6 ± 23.9	37.2 ± 19.1	30.2 ± 19.6
			(P<0.01)	(P=0.004)
	Hemoglobin, g/dL	8.4 ± 1.4	8.6 ± 1.6	8.5 ± 1.8
			(P=0.568)	(P=0.386)
	Serum creatinine, mg/dL	0.7 ± 0.2	0.8 ± 0.3	0.8 ± 0.4
			(P=0.09)	(P=0.124)
	Cystatin-C, mg/L	0.9 ± 0.4	1.1 ± 0.4	1.2 ± 0.6
			(P=0.001)	(P<0.001)
	Prot 24h, mg/24h	440.4 ± 578.9	477.2 ± 450.11	469.6 ± 466.7
			(P=0.065)	(P=0.073)
	NM, value not measured
	HIC, hepatic iron concentration
	LVEF, left ventricular ejection fraction

Abstract

Iron overload is present in several cases of double heterozygous sickle-cell/beta-thalassemia (HbS/thal). Deferasirox is an orally administered iron chelator which is effective on iron overloaded patients with transfusion-dependent anemia. The aim of this study was to investigate the efficacy and safety of deferasirox on HbS/thal patients with iron overload.

We evaluated 31 adult patients with HbS/thal (14M/17F; median age 41 years) who had serum ferritin levels >1000 ng/mL and who were sporadically transfused. Total iron burden was monitored by measuring serum ferritin levels before and monthly after starting deferasirox, while liver iron concentration and cardiac iron burden were measured by magnetic resonance imaging (MRI) T2 and T2* parameters at baseline and 12 months after deferasirox treatment.

Deferasirox managed to reduce the mean serum ferritin levels after 12 months of treatment from 1989±923 to 1008±776 ng/mL (P<0.001). This reduction was accompanied by a significant improvement on MRI T2* of the liver (from 3.9±3.2 to 5.8±3.1 ms; P<0.01) and by a comparable improvement of biochemical parameters of liver function. Mild nausea and diarrhoea of grade 1/2 were reported in 25% of patients within the first month of treatment, but did not reoccur. These data indicate that deferasirox provided effective control of iron levels (mainly of the liver) in minimally transfused patients with HbS/thal, without significant adverse events, at similar doses to those studied widely for the treatment of patients with thalassemia syndromes.

Legend to Figure

Deferasirox produced a dramatic reduction of serum ferritin levels after 12 months of treatment (see also Table 1)


Tilemahos J Kalvouridis
email: tkalvouridis@gmail.com

Stationary solutions of a small gyrostat in the Newtonian field of two bodies with equal masses

Keywords: Gyrostat, Rigid body dynamics, Stationary solutions, Stability

The paper deals with the dynamics of a small gyrostat satellite acted upon by the Newtonian forces of two big bodies of equal masses which rotate around their center of mass. The gyrostat's equations of motion are derived and classes of its stationary solutions, as well as their stability are studied.

Introduction

Problems with gyrostats have been studied in the past, particularly during the last decades ([1-10, 12-14], etc.). Recently, Vera and Vigueras [START_REF] Vera | Libration points of a gyrostat in a generalized restricted three-body problem: linear stability[END_REF][START_REF] Vera | Hamiltonian dynamics of a gyrostat in the n-body problem: relative equilibria[END_REF] and Vera [START_REF] Vera | Eulerian equilibria of a triaxial gyrostat in the three-body problem: rotational Poisson dynamics in Eulerian equilibria[END_REF] have renewed the interest on this issue with important results. These problems are extremely interesting, since most of the artificial satellites and spacecrafts have mobile parts such as rotating antennas, photovoltaic panels, or mechanical articulated arms, etc., and can therefore be presumed as gyrostatic bodies. On the other hand, a special case of the old and famous restricted three-body problem is the Copenhagen configuration where the two dominant masses are equal. This case has lately gained new scientific interest after the discovery of many exo-solar planetary systems, the majority of which consists of two members. The recent estimations about the partners in such a system show that, in many cases, their masses can roughly be considered as equal. Here, we consider a gyrostat S with infinitesimal mass moving in the Newtonian field of two big spherical and homogeneous bodies P 1 and P 2 of equal masses m that rotate about their center of mass with constant angular velocity ω. The much smaller body S of mass m 0 , does not affect the motion of the primaries and has a gyrostatic structure (Fig. 1). As we know, a gyrostat consists of two parts: a rigid part called the platform or carrier and n other parts, the rotors, which are connected to the platform. Each rotor is spinning-independently about an axis fixed on the platform and its motion does not modify the mass distribution of the gyrostat. The platform may rotate about an inertial reference frame, so that a gyrostat is generally characterized by n + 1 angular velocities, n of which, are the angular velocities of the spinning rotors relative to the platform. Sometimes, in the relevant literature, a gyrostat is also referred to as a 'dual-spin' body. Hereunder, we shall examine gyrostats with only one rotor. First, we extract the equations of motion of the gyrostat and then we investigate the equilibrium states (translational and rotational) and their linear stability assuming that the small body is inertially axisymmetric. The obtained results are illustrated with some figures.

Fig. 1 The configuration of the system. The framed picture shows the gyrostatic structure of the infinitesimal body S and the body-fixed axes. GN is the line of nodes and Oxyz is the synodic coordinate system

Equations of motion

We use a synodic system Oxyz (unit vectors i, j, k) rigidly attached to the system of the two primaries P 1 and P 2 and centered at their mass center O. Its xaxis coincides with the line joining the two primaries, while the y-axis is perpendicular to it on the orbital plane of the massive bodies. We assume that the synodic system rotates with a constant angular velocity ω. We also use two auxiliary reference frames centered at the mass center G of S; the body-fixed GXY Z (unit vectors I , J , K) the axes of which coincide with the gyrostat's principal central axes of inertia and Gx y z , the axes of which are parallel to those of the synodic system (Fig. 1).

In our calculations, we consider the triad of Euler angles (ψ, θ, φ) that is defined from the rotating sequence 3-1-3 (Fig. 1). Hence, the components of the angular velocity Ω of the gyrostat with respect to an inertial frame, can be expressed as follows:

Ω = ( ψ + ω) sin φ sin θ + θ cos φ, ( ψ + ω) cos φ sin θ -θ sin φ, ( ψ + ω) cos θ + φ
The kinetic energy of the translational-rotational motion of the gyrostat is given with the relation:

T = 1 2 m 0 ( ẋ -ωy) 2 + ( ẏ + ωx) 2 + ż2 + 1 2 (I 1 -I 2 ) sin θ sin 2φ θ( ψ + ω) + 1 2 I 1 sin 2 φ + I 2 cos 2 φ sin 2 θ + I 3 cos 2 θ ( ψ + ω) 2 + 1 2 I 1 cos 2 φ + I 2 sin 2 φ θ 2 + 1 2 I 3 φ2 + I 3 cos θ φ( ψ + ω) + (h a cos φ -h b sin φ) θ + c φ + (h a sin φ + h b cos φ) sin θ + c cos θ ( ψ + ω) + T r
where m 0 is the mass of the gyrostat, I i , i = 1, 2, 3 are its principal central moments of inertia, h a , h b , h c are the components of the rotor's internal moment of momentum and T r is its kinetic energy relative to the platform. If we assume that the rotor rotates with relation to the carrier with constant angular velocity, then the quantity T r is constant, and consequently it does not appear in the final form of the equations of motion.

The potential energy after MacCullagh's formula [START_REF] Roy | Orbital Motion[END_REF] can be expressed in terms of Eulerian angles [START_REF] Michalakis | The equilibrium configurations of the restricted problem of 2 + 2 triaxial rigid bodies[END_REF]

V = -Gmm 0 1 r 1 + 1 r 2 + Gm 2 i=1 1 r 5 i (x -x i ) 2 f + y 2 g + z 2 u + (x -x i )yh + (x -x i )zw + yzv
where r 1 , r 2 , are the distances of the gyrostat S from the primaries,

r 2 1 = (x -x 1 ) 2 + y 2 + z 2 r 2 2 = (x -x 2 ) 2 + y 2 + z 2 and f = - 3p 2
p = I 1 + I 2 -2I 3 2 and s = I 1 -I 2
The equations of motion are of two kinds: those that express the translational motion of the body and those that describe its rotational motion. Both systems of equations are coupled. More specifically:

I. Translational motion The translational motion is described by means of the equations

ẍ -2ω ẏ = ω 2 x -Gm 2 i=1 (x -x i ) r 3 i -G m m 0 2 i=1 ∂ ∂x 1 r 5 i (x -x i ) 2 f + y 2 g + z 2 u + (x -x i )yh + (x -x i )zw + yzv (1) ÿ + 2ω ẋ = ω 2 y -Gm 2 i=1 y r 3 i -G m m 0 2 i=1 ∂ ∂y 1 r 5 i (x -x i ) 2 f + y 2 g + z 2 u + (x -x i )yh + (x -x i )zw + yzv (2) z = -Gm 2 i=1 z r 3 i -G m m 0 2 i=1 ∂ ∂z 1 r 5 i (x -x i ) 2 f + y 2 g + z 2 u + (x -x i )yh + (x -x i )zw + yzv (3)
II. Rotational motion The rotational motion is described by the equations

I 1 sin 2 φ + I 2 cos 2 φ sin 2 θ + I 3 cos 2 θ ψ + 1 2 (I 1 -I 2 ) sin 2φ sin θ θ + I 3 cos θ φ + (I 1 -I 2 ) sin 2φ sin 2 θ φ + I 1 sin 2 φ + I 2 cos 2 φ sin 2θ θ -I 3 sin 2θ θ ( ψ + ω) + (I 1 -I 2 ) cos 2φ sin θ θ φ + 1 2 (I 1 -I 2 ) sin 2φ cos θ θ 2 -I 3 sin θ θ φ + (h a sin φ + h b cos φ) cos θ -h c sin θ θ + (h a cos φ -h b sin φ) sin θ φ = -Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂ψ + y 2 ∂g ∂ψ + z 2 ∂u ∂ψ + (x -x i )y ∂h ∂ψ + (x -x i )z ∂w ∂ψ + yz ∂v ∂ψ (4) 1 2 (I 1 -I 2 ) sin 2φ sin θ ψ + I 1 cos 2 φ + I 2 sin 2 φ θ -(I 1 -I 2 ) sin 2φ θ φ + (I 1 -I 2 ) cos 2φ sin θ φ( ψ + ω) - 1 2 I 1 sin 2 φ + I 2 cos 2 φ -I 3 × sin 2θ( ψ + ω) 2 + I 3 sin θ φ( ψ + ω) -(h a sin φ + h b cos φ) φ -(h a sin φ + h b cos φ) cos θ -h c sin θ × ( ψ + ω) = -Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂θ + y 2 ∂g ∂θ + z 2 ∂u ∂θ + (x -x i )y ∂h ∂θ + (x -x i )z ∂w ∂θ + yz ∂v ∂θ ( 5
)
I 3 cos θ ψ + I 3 φ -I 3 sin θ θ( ψ + ω) - 1 2 (I 1 -I 2 ) sin 2φ sin 2 θ( ψ + ω) 2 + 1 2 (I 1 -I 2 ) sin 2φ θ 2 -(I 1 -I 2 ) cos 2φ × sin θ θ( ψ + ω) -(h a cos φ -h b sin φ) × sin θ( ψ + ω) + (h a sin φ + h b cos φ) θ = -Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂φ + y 2 ∂g ∂φ + z 2 ∂u ∂φ + (x -x i )y ∂h ∂φ + (x -x i )z ∂w ∂φ + yz ∂v ∂φ (6)
As we have mentioned before, (1) to (6) are coupled.

Stationary solutions and stability

Stationary solutions of the gyrostat exist under the conditions ẋ = ẏ = ż = ẍ = ÿ = z = θ = θ = φ = φ = ψ = ψ = 0. For practical reasons we shall confine our study to those orientations which are determined by all possible combinations of the angles

ψ = λ ψ π 2 with λ ψ = 0, 1, 2, 3 θ = λ θ π 2 with λ θ = 0, 1, 2 and φ = λ φ π 2 with λ φ = 0, 1, 2, 3
For these values, it holds that

h = v = w = 0
Then (3) takes the form

-Gmm 0 z 1 r 3 1 + 1 r 3 2 -Gm 2 i=1 -5z(x -x i ) 2 f r 7 i - 5zy 2 g r 7 i - 5z 3 u r 7 i + 2z
r 5 i u = 0 or -Gmm 0 z 1 r 3 1 + 1 r 3 2 -Gm 2 i=1 z r 7 i -5(x -x i ) 2 f -5y 2 g -5z 2 u + 2r 2 i u = 0 and finally z -m 0 1 r 3 1 + 1 r 3 2 - 2 i=1 1 r 7 i -5(x -x i ) 2 f -5y 2 g -5z 2 u + 2r 2 i u = 0 ( 7 )
This relation proves the existence of solutions on xy plane (z = 0). Hereunder, we shall confine our study to the planar equilibrium positions. Then

-Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂ψ + y 2 ∂g ∂ψ + (x -x i )y ∂h ∂ψ = 0 ( 8 ) -Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂θ + y 2 ∂g ∂θ + (x -x i )y ∂h ∂θ + (h a sin φ + h b cos φ) cos θ -h c sin θ ω = 0 ( 9 ) -Gm 2 i=1 1 r 5 i (x -x i ) 2 ∂f ∂φ + y 2 ∂g ∂φ + (x -x i )y ∂h ∂φ + (h a cos φ -h b sin φ) sin θ ω = 0 ( 10
) Table 1 Case 1 (I 1 = I 2 ) I 1 = I 2 Group h a = 0, h b = 0, h c = 0 Group h a = 0, h b = 0, h c = 0 ψ θ φ ψ θ φ I 0 0 π 2 V 0 0 0 0 0 3π 2 0 0 π π 0 π 2 π 0 0 π 0 3π 2 π 0 π II π 2 0 π 2 VI π 2 0 0 π 2 0 3π 2 π 2 0 π 3π 2 0 π 2 3π 2 0 0 3π 2 0 3π 2 3π 2 0 π III 0 π π 2 VII π 2 π 0 0 π 3π 2 π 2 π π π π π 2 3π 2 π 0 π π 3π 2 3π 2 π π IV π 2 π π 2 VIII 0 π 0 π 2 π 3π 2 0 π π 3π 2 π π 2 π π 0 3π 2 π 3π 2 π π π
It is easily proved that when the Euler angles are integer multiples of π/2, then

∂f ∂ψ = ∂f ∂θ = ∂f ∂φ = ∂g ∂ψ = ∂g ∂θ = ∂g ∂φ = ∂h ∂θ = 0
and in this case,

-Gm 2 i=1 1 r 5 i (x -x i )y ∂h ∂ψ = 0 (11) (h a sin φ + h b cos φ) cos θ -h c sin θ ω = 0 (12) -Gm 2 i=1 1 r 5 i (x -x i )y ∂h ∂φ + (h a cos φ -h b sin φ) sin θ ω = 0 ( 13
)
The above [START_REF] Roy | Orbital Motion[END_REF] to [START_REF] Tong | Chaotic motion of an asymmetric gyrostat in the gravitational field[END_REF] are verified for any triad (ψ, θ, φ) of Euler angles that are integer multiples of π/2 provided that the conditions written at the top of each table, are valid. We observe that each condition, I 1 = I 2 , I 1 = I 3 , or I 2 = I 3 implies f = g. More specifically, we have the following cases:

Table 2 Case 2 (I 1 = I 3 ) I 1 = I 3 h a = 0, h b = 0, h c = 0 ψ θ φ 0 π 2 0 0 π 2 π π 2 π 2 0 π 2 π 2 π π π 2 0 π π 2 π 3π 2 π 2 0 3π 2 π 2 π Case 1. I 1 = I 2 . Then f = g = 1 2 (I 2 -I 3
). Equations [START_REF] Roy | Orbital Motion[END_REF] to [START_REF] Tong | Chaotic motion of an asymmetric gyrostat in the gravitational field[END_REF] are verified for the combinations of the Euler angles which appear in Table 1. Case 2. I 1 = I 3 . Then f = g = 1 2 (I 1 -I 2 ). Equations ( 11) to [START_REF] Tong | Chaotic motion of an asymmetric gyrostat in the gravitational field[END_REF] are verified for the combinations of the Euler angles which appear in Table 2. Case 3. I 2 = I 3 . Then f = g = 1 2 (I 3 -I 1 ) and ( 11)-( 13) are verified for the combinations of the Euler angles which appear in Table 3.

Table 3 Case 3 (I 2 = I 3 ) I 2 = I 3 h a = 0, h b = 0, h c = 0 ψ θ φ 0 π 2 π 2 0 π 2 3π 2 π 2 π 2 π 2 π 2 π 2 3π 2 π π 2 π 2 π π 2 3π 2 3π 2 π 2 π 2 3π 2 π 2 3π 2
In view of the above values of f and g, ( 1) and ( 2) take the form

ω 2 x -Gm 2 i=1 (x -x i ) r 3 i -3G m m 0 f 2 i=1 (x -x i ) r 5 i = 0 ( 14
)
ω 2 y -Gmy 2 i=1 1 r 3 i -3G m m 0 yf 2 i=1 y r 5 i = 0 ( 15
)
We shall transform the above equations to dimensionless ones, in order to show their relation to those of a point-like mass. Thus, we first replace the principal central moments of inertia I i , i = 1, 2, 3 with their corresponding radii of gyration

ρ 2 i = I i m 0 , i = 1, 2, 3. Then we have f = 1 2 m 0 ρ 2 2 -ρ 2 3 for case 1 f = 1 2 m 0 ρ 2 1 -ρ 2 2
for case 2 and

f = 1 2 m 0 ρ 2 3 -ρ 2 1 for case 3 Provided that a = |P 1 P 2 |, ω 2 = Gm a 3
and by introducing the transformations

x = ax * , y= ay * , r i = ar * i r = ar * and f = a 2 m 0 (f * )
Fig. 2 The considered configuration with an inertially axisymmetric gyrostat (I 1 = I 2 ) when h b = 0. The framed pictures show three different orientations of attitude equilibrium obtained from Table 1 into ( 14) and ( 15), we obtain the dimensionless relations

x -

2 i=1 (x -x i ) r 3 i + 3 2 f 2 i=1 (x -x i ) r 5 i = 0 (16) y - 2 i=1 y r 3 i + 3 2 f 2 i=1 y r 5 i = 0 (17)
from which we have omitted the asterisks. The orbital equilibria of m 0 with the constant orientations of the body-fixed axes are the solutions of the system of ( 16) and [START_REF] Vera | Eulerian equilibria of a triaxial gyrostat in the three-body problem: rotational Poisson dynamics in Eulerian equilibria[END_REF]. It is evident that when S has a spherical symmetry and the rotor is motionless in relation to the platform, then it will be

I 1 = I 2 = I 3 , h a = h b = h c = 0 f = g = h = u = v = w = 0
and the problem reduces to the classic Copenhagen case where S is a material point. Figure 2 shows some stationary solutions of an inertially axisymmetric gyrostat S with I 1 = I 2 . The gyrostat rests on an orbital equilibrium position. The framed pictures show three different orientations of attitude equilibrium. The correspondence of the displayed orientations with Euler angles is obtained easily from Table 1.

Let us assume that L e (x e , y e , 0, ψ e , θ e , φ e ) is an equilibrium state of the gyrostat. If the conditions that correspond to this position are slightly disturbed, the variational equations of the planar translational motion will take the form:

δ ẍ -2ωδ ẏ = q=x,y,ψ,θ,φ ∂E 1 ∂q L e δq δ ÿ + 2ωδ ẋ = q=x,y,ψ,θ,φ ∂E 2
∂q L e δq E 1 and E 2 are the right parts of ( 1) and ( 2) while their time derivatives are calculated on the position L e .

Similarly, the equations of the rotational motion become In the last three equations, we denote with E i , i = 3, 4, 5 the right parts of (4), [START_REF] El-Gohary | On the stability of relative programmed motion of satellite-gyrostat[END_REF], and (6), while their time derivatives are calculated on the equilibrium point L e . By considering orientations ψ e , θ e , φ e that are integer multiples of π/2, it is easily proven that

∂E i ∂q L e = 0, i = 1, 2, q = ψ, θ, φ and ∂E j ∂Q L e = 0, i = 1, 2, Q = x, y
Then the variational equations can be written as We note that the variational equations that are related to the translation of the gyrostat, as well as those related to its rotation are uncoupled between them. Studying the stability of the gyrostat's rotational equilibria-for the various combinations of Euler angles accompanied by the corresponding conditions between I 1 , I 2 , I 3 and h a , h b , h c as classified on Tables 1,2, and 3-it can be easily proven that for the cases considered it will always be δψ(t) = δ ψe t + δψ e (23)

δ ẍ -2ωδ ẏ = ∂E 1 ∂x L e δx
If, for example, we consider the triad ψ = 0, θ = π/2, φ = π/2of Table 3, then (20) becomes

I 1 δ ψ = 0
the solution of which has form (23). From all these, we come to the conclusion that the deviation of the disturbed solution is unbounded and consequently all the above equilibrium states of the gyrostat are unsteady.

Conclusions

The principal conclusion derived from the above analysis is that the stationary solutions of the small gyrostat which correspond to the integer multiple of π/2 of the Eulerian angles are all unstable. Closing this article, we would like to emphasize anew the value of this model, especially as regards its ability to approach physical systems as described in detail in the introduction. In addition, stability appears as a necessary part of the equilibrium solutions. Especially in Gyro-dynamics the conclusions are quite important regardless of the actual character of equilibria (stable or unstable), because of their involvement to practical applications such as the transition of a satellite from an equilibrium position to another, the re-orientation of a vehicle, etc.

I 1 2 (I 1 -I 2 ) 2 (I 1 -I 2 ) 2 (I 1 - 2 I 1 2 (I 1 - 2 (I 1 -

 121221221212121 sin 2 φ e + I 2 cos 2 φ e sin 2 θ e + I 3 cos 2 θ e δ ψ + 1 sin 2φ e sin θ e δ θ + I 3 cos θ e δ φ + I 1 sin 2 φ e + I 2 cos 2 φ e -I 3 ω sin 2θ e + (h a sin φ e + h b cos φ e ) cos θ eh c sin θ e δ θ + (I 1 -I 2 )ω sin 2 θ e sin 2φ e + (h a cos φ e h b sin φ e ) sin θ e δ φ = q=x,y,ψ,θ,φ ∂E 3 ∂q L e δq 1 sin 2φ e sin θ e δ ψ + I 1 cos 2 ϕ e + I 2 sin 2 ϕ e δ θ -I 1 sin 2 φ e + I 2 cos 2 φ e -I 3 ω sin 2θ e + (h a sin φ e + h b cos φ e ) cos θ e h c sin θ e δ ψ + (I 1 -I 2 )ω cos 2φ e sin θ e + I 3 ω sin θ e -(h a sin φ e + h b cos φ e ) δ φ -I 1 sin 2 φ e + I 2 cos 2 φ e -I 3 ω 2 cos 2θ e -(h a sin φ e + h b cos φ e ) sin θ e + h c cos θ e ω × δθ -1 I 2 )ω 2 sin 2θ e sin 2φ e + (h a cos φ eh b sin φ e )ω cos θ e δφ -1 sin 2 φ e + I 2 cos 2 φ e -I 3 ω 2 sin 2θ e -(h a sin φ e + h b cos φ e ) cos θ eh c sin θ e ω = q=x,y,ψ,θ,φ ∂E 4 ∂q L e δq I 3 cos θ e δ ψ + I 3 δ φ -(I 1 -I 2 )ω sin 2φ e sin 2 θ e + (h a cos φ eh b sin φ e ) sin θ e δ ψ -(I 1 -I 2 )ω cos 2ϕ e sin θ e + I 3 ω sin θ e -(h a sin φ e + h b cos φ e ) δ θ -1 I 2 )ω 2 sin 2θ e sin 2φ e + (h a cos φ eh b sin φ e )ω cos θ e δθ -(I 1 -I 2 )ω 2 cos 2φ e sin 2 θ e -(h a sin φ e + h b cos φ e )ω sin θ e δφ -1 I 2 )ω 2 sin 2φ e sin 2 θ e + (h a cos φ e h b sin φ e )ω sin θ e = q=x,y,ψ,θ,φ ∂E 5 ∂q L e δq

I 3

 3 cos θ e δ ψ + I 3 δ φ -(h a cos φ eh b sin φ e ) sin θ e δ ψ -(I 1 -I 2 )ω cos 2φ e sin θ e + I 3 ω sin θ e -(h a sin φ e + h b cos φ e ) δ θ -(h a cos φ eh b sin φ e )ω cos θ e δθ -(I 1 -I 2 )ω 2 cos 2φ e sin 2 θ e -(h a sin φ e + h b cos φ e )ω sin θ e δφ -(h a cos φ eh b sin φ e )ω sin θ e

  I 1 sin 2 φ e + I 2 cos 2 φ e sin 2 θ e + I 3 cos 2 θ e δ ψ + I 3 cos θ e δ φ + (h a sin φ e + h b cos φ e ) cos θ e h c sin θ e δ θ + (h a cos φ eh b sin φ e ) sin θ e δ φ I 1 cos 2 φ e + I 2 sin 2 φ e δ θ + (h a sin φ e + h b cos φ e ) cos θ eh c sin θ e δ ψ + (I 1 -I 2 )ω cos 2φ e sin θ e + I 3 ω sin θ e -(h a sin φ e + h b cos φ e ) δ φ -I 1 sin 2 φ e + I 2 cos 2 φ e -I 3 ω 2 cos 2θ e -(h a sin φ e + h b cos φ e ) sin θ

	=	q=ψ,θ,φ	∂E 3 ∂q L e	δq	(20)
	=	q=ψ,θ,φ	∂E 4 ∂q L e	δq	(21)
							+	∂E 1 ∂y L e	δy	(18)
					δ ÿ + 2ωδ ẋ =	∂E 2 ∂x L e	δx +	∂E 2 ∂y L e	δy	(19)

e + h c cos θ e ω δθ -(h a cos φ e h b sin φ e )ω cos θ e δφ -(h a sin φ e + h b cos φ e ) cos θ eh c sin θ e ω