Peter T Kim
email: pkim@uoguelph.ca.

Donald P St

⋆⋆ Richards
email: richards@stat.psu.edu.

Diffusion Tensor Imaging and Deconvolution on Spaces of Positive Definite Symmetric Matrices

Keywords: Deconvolution, diffusion tensor imaging, Harish-Chandra c-function, Helgason-Fourier transform, inverse problems, Laplace-Beltrami operator, magnetic resonance imaging, Sobolev spaces, Wishart distribution

Diffusion tensor imaging can be studied as a deconvolution density estimation problem on the space of positive definite symmetric matrices. We develop a nonparametric estimator for the common density function of a random sample of positive definite matrices. Our estimator is based on the Helgason-Fourier transform and its inversion, the natural tools for analysis of compositions of random positive definite matrices. Under smoothness conditions on the density of the intrinsic error in the random sample, we derive bounds on the rates of convergence of our nonparametric estimator to the true density.

Introduction

The appearance in medical imaging of sequences of random positive definite matrices has become commonplace due to developments in diffusion tensor imaging (DTI), a method of imaging based upon the observation that water molecules in biological tissue are always in motion. For the purposes of mathematical modeling, it is generally assumed that the diffusion of water molecules at any given location in biological tissue follows a Brownian motion. A diffusion tensor image then is represented by the 3 × 3 positive-definite covariance matrix of the local diffusion process at the given location (Fletcher and Joshi [START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF][START_REF] Fletcher | Riemannian geometry for the statistical analysis of diffusion tensor data[END_REF]. DTI seeks to detect the diffusion of water protons between and within distinct tissue cells, and to derive estimates of the dominant orientation and direction of the Brownian motion (Le Bihan [START_REF] Le Bihan | Diffusion and Perfusion Magnetic Resonance Imaging: Applications to Functional MRI[END_REF], Hasan, et al. [START_REF] Hasan | Analytical computation of the eigenvalues and eigenvectors in DT-MRI[END_REF]).

In DTI brain imaging, the diffusion of water molecules within and between voxels, the three-dimensional volume elements that constitute an image, reveal both the orientation of fibers comprising white-matter tracts in the brain and the coherence of fibers, the extent to which fibers are aligned together. DTI may be the only non-invasive, in vivo procedure which enables the study of deep brain white-matter fibers. Consequently, DTI has been found to be highly promising for comparing the human brain in normal states with abnormal states caused by strokes, epileptic seizures, tumors, white-matter abnormalities, multiple sclerosis lesions, HIV-infection, traumatic brain injuries, aging, Alzheimer's disease, alcoholism, and developmental disorders; and there are potential applications to psychiatric conditions including schizophrenia, autism, cognitive and learning disabilities (Neumann-Haefelin, et al. [START_REF] Neumann-Haefelin | New magnetic resonance imaging methods for cerebrovascular disease[END_REF], Rosenbloom, et al. [START_REF] Rosenbloom | Using magnetic resonance imaging and diffusion tensor imaging to assess brain damage in alcoholics[END_REF], Pomara, et al. [START_REF] Pomara | White-matter abnormalities in HIV-1 infection: A diffusion tensor imaging study[END_REF], Matthews and Arnold [START_REF] Matthews | Magnetic resonance imaging of multiple sclerosis: new insights linking pathology to clinical evolution[END_REF]). In addition, DTI has been applied in research on the pathology of organ and tissue types such as the human breast, kidney, lingual, cardiac, skeletal muscles, and spinal cord (Damon,et al. [2]).

It is well-known that magnetic resonance imaging, from which diffusion tensors are derived, is endowed inhererently with random noise. Hence, DTI data also contain noise (Basu,et al. [1]), and it is natural that statistical inferential issues arise in the analysis of DTI data (Koltchinskii,et al. [8]; Schwartzman, et al. [START_REF] Schwartzman | Random Ellipsoids and False Discovery Rates: Statistics for Diffusion Tensor Imaging Data[END_REF][START_REF] Schwartzman | Cross-subject comparison of principal diffusion direction maps[END_REF]; Zhu, et al. [START_REF] Zhu | Statistical analysis of diffusion tensors in diffusion-weighted magnetic image resonance data (with discussion)[END_REF]).

In this paper, we study the problem of estimating the probability density function of a population of positive definite matrices based on a random sample from that population. An instance in which this problem arises may be obtained from Schwartzman [START_REF] Schwartzman | Random Ellipsoids and False Discovery Rates: Statistics for Diffusion Tensor Imaging Data[END_REF] who studied the two-sample comparison of twelve children divided into two groups according to reading ability, where the issue is to compare physical characteristics of brain tissue of the two groups on the basis of DTI images. In addition to comparing the population parameters, it is natural to seek estimators of the underlying density functions, and then it will be important to estimate the rates of convergence of the density estimators.

The deconvolution density estimation problem has been widely studied on Euclidean spaces. this classical setting, the commutative nature of the underlying mathematical operations renders the problem amenable to classical mathematical methods. The deconvolution problem has also been studied on certain compact manifolds. In that setting, the problem is solvable using well-known generalizations of classical Fourier analysis. By contrast, P m , the space of m × m positive definite matrices, is a noncompact Riemannian manifold and has an intrinsic non-commutative nature, and it is natural to expect that the deconvolution problem will be more difficult in that setting. To the best of our knowledge, no results yet are available for the general deconvolution problem on P m .

To date, the primary statistical emphases regarding mathematical methods on P m are motivated by properties of the Wishart distribution, see for example Muirhead [START_REF] Muirhead | Aspects of Multivariate Statistics[END_REF]. To solve the general deconvolution problem on P m , more advanced mathematical methods are required. While much is known in the mathematical literature about the necessary methods (Helgason [6], Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF]), virtually nothing about those methods has appeared in the statistical literature. Using these new methods, we develop nonparametric methods for solving the deconvolution problem on P m . In summary, section 2 provides notation, and introduces the Helgason-Fourier transform and its inversion formula, involving the Harish-Chandra c-function; we have provided the necessary details of these concepts so as to make the paper fully accessible to readers who are new to this area. We formalize in section 3 the statistical procedure in terms of measurement errors on P m and present the main results on general deconvolution density estimation on P m and we provide the explicit details in the case of the Wishart distribution. All proofs can be found in Kim and Richards [START_REF] Kim | Deconvolution density estimation on spaces of positive definite symmetric matrices[END_REF].

Preliminaries

Throughout the paper, we denote by G the general linear group GL(m, R), of m×m nonsingular real matrices and by K the group, O(m), of m×m orthogonal matrices. The group G acts transitively on P m , the space of m × m positive definite matrices, by the action

G × P m → P m , (g, x) → g ′ xg,
g ∈ G, x ∈ P m , where g ′ denotes the transpose of g. Under this action, the isotropy group of the identity in G is K, hence the homogeneous space K\G can be identified with P m by the "natural" mapping

K\G → P m , Kg → g ′ g. (1)
In distinguishing between left and right cosets, we place the quotient operation on the left and right of the group, respectively. A random matrix X ∈ P m is said to be K-invariant if

X L = k ′ Xk for all k ∈ K, where " L =" denotes equality in distribution. A function f on P m is called K-invariant if f (k ′ xk) = f (x)
for all k ∈ K, x ∈ P m ; we will indicate that f is K-invariant by writing its domain as P m /K, with a similar notation for K-invariant positive definite random matrices.

By means of the relationship (1) between K\G and P m , we identify Kinvariant functions on P m with K-biinvariant functions on G, i.e., functions

f : G → C which satisfy f (g) = f (k 1 gk 2 ) for all k 1 , k 2 ∈ K and g ∈ G.
In particular, P m /K ≃ K\G/K where "≃" denotes diffeomorphic equivalence.

Consider random matrices X, ε ∈ P m with corresponding group elements X, ε ∈ G, respectively. By the natural map (1), K X → X, equivalently X ′ X = X and, similarly,

K ε → ε. Then X ε ∈ G is mapped via (1) to X ε → ( X ε) ′ ( X ε) ≡ ε ′ X ′ X ε = ε ′ X ε ∈ P m . ( 2
) If ε is K-biinvariant, i.e., ε L = k 1 εk 2 for all k 1 , k 2 ∈ K then ε ′ ε = ε = k ′ Λk where k ∈ K and Λ ∈ P m is the diagonal matrix of eigenvalues of ε. Consequently, we define ε 1/2 = k ′ Λ 1/2 k
and note that, by K-biinvariance of ε, the relationship (2) on the group G corresponds in distribution to ε 1/2 Xε 1/2 ∈ P m , where X ∈ P m and ε ∈ P m /K. Bearing this in mind, we formally make the following definition. Suppose that X ∈ P m and ε ∈ P m /K are random matrices. Then the composition of X and ε is

X • ε = ε 1/2 Xε 1/2 (3)
where ε 1/2 is the positive definite square root of ε.

There is an alternative approach leading to compositions (3). We begin by noting that P m is a Riemannian manifold, hence each pair X, Y ∈ P m defines a unique geodesic path. Since P m is also a homogeneous space, viz., P m = K\G, then there exists a unique V ∈ G such that

Y = V XV ′ . ( 4
)
This model can be viewed as a multiplicative analog of classical regression analysis, with Y serving as the dependent variable, X as the independent variable, and V as the error variable. We can also study this relationship between X and Y using the logarithm map on P m . To that end, each observation on Y may be transformed into an observation on y where exp(y) = Y with y ∈ T (P m ), the tangent space of P m . Similarly, each measurement on X may be transformed into a measurement on x ∈ T (P m ) where exp(x) = X. By postulating the existence of a "small" error v ∈ T (P m ) such that y = x + v and, by exponentiating this linear "regression" relationship between y and x, we are led naturally to (4); see Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], section 4.1-4.2 for details on the exponentiation of such linear relationships on T (P m ).

If we assume that the measurement error V is isotropic or has no preferred orientation, i.e., V has a K-invariant distribution, then the conclusion is the model ( 4) in which V ∈ G/K. As noted in [START_REF] Rosenbloom | Using magnetic resonance imaging and diffusion tensor imaging to assess brain damage in alcoholics[END_REF], in the context of DTI, the assumption that water molecules diffuse istropically is appropriate for regions such as the ventricles, which are large fluid-filled spaces deep in the brain. On the other hand, water molecules located in white-matter fiber are constrained by the axon sheath; this forces greater movement along the longitudinal axes of fibers than across the axes, and then diffusion may be isotropic only at sufficiently small scales.

Returning to (4), we apply polar coordinates on G, viz.,

V = kε 1/2 where k ∈ K and ε ∈ P m . Since V is K-invariant then V L = k ′ V = k ′ kε 1/2 = ε 1/2 ; hence (4) reduces to Y L = ε 1/2 Xε 1/2 ≡ X •ε, which agrees with (3) since ε = V ′ V is K-invariant.
Thus, the problem is to estimate nonparametrically the density of X based on a random sample from Y = X • ε where the error matrix ε is isotropic.

Let C ∞ c (P m ) denote the collection of complex-valued, infinitely differentiable, compactly supported functions f on P m . For w ∈ P m and j = 1, . . . , m, denote by |w j | the principal minor of order j of w. For s = (s 1 , . . . , s m ) ∈ C m , the power function p s : P m → C is defined by

p s (w) = m j=1 |w j | sj , w ∈ P m .
Let dk denote the Haar measure on K, normalized to have total volume equal to one. Then

h s (w) = K p s (k ′ wk) dk (5)
w ∈ P m , s ∈ C m , is a zonal spherical function on P m . It is well-known that the spherical functions play a fundamental role in harmonic analysis on symmetric spaces; see Helgason [START_REF] Helgason | Groups and Geometric Analysis[END_REF]. In particular, ( 5) is a special case of Harish-Chandra's formula for the general spherical functions. If s 1 , . . . , s m are nonnegative integers then, except for a constant factor, ( 8) is an integral formula for the zonal polynomials which arise often in aspects of multivariate distribution theory; see Muirhead [START_REF] Muirhead | Aspects of Multivariate Statistics[END_REF], pp. 231-232. Let w = (w ij ) ∈ P m . Then, up to a constant factor, the unique G-invariant

measure on P m is d * w = |w| -(m+1)/2 1≤i≤j≤m dw ij , (6)
where |w| is the determinant of w.

For s = (s 1 , . . . , s m ) ∈ C m and k ∈ K, the Helgason-Fourier transform

([17], p. 87) of a function f ∈ C ∞ c (P m ) is Hf (s, k) = Pm f (w) p s (k ′ wk) d * w, (7)
where p s (k ′ wk) denotes complex conjugation and d * w is the G-invariant measure (6).

For the case in which f ∈ C ∞ c (P m /K), we make the change of variables w → k ′ 1 wk 1 in (7), k 1 ∈ K, and integrate with respect to the Haar measure dk 1 . Applying the invariance of f and the formula (5) for the zonal spherical function, we deduce that Hf (s, k) does not depend on k; specifically, Hf (s, k) = f (s) where

f (s) = Pm f (w) h s (w) d * w, (8)
is the zonal spherical transform of f . Let A = {diag(a 1 , . . . , a m ) : a j > 0, j = 1, . . . , m}, denote the group of diagonal positive definite matrices in G, and N = {n = (n ij ) ∈ G : n ij = 0, 1 ≤ j < i ≤ m; n jj = 1, j = 1, . . . , m} be the subgroup of G consisting of uppertriangular matrices with all diagonal entries equal to 1. It is well-known (see Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], p. 20) that each g ∈ G can be decomposed uniquely as g = kan where k ∈ K, a ∈ A, and n ∈ N ; this result is the Iwasawa decomposition, and (k, a, n) are called the Iwasawa coordinates of g.

c m (s) = 1≤i<j≤m-1 B 1 2 , s i + • • • + s j + 1 2 (j -i + 1) B 1 2 , 1 2 (j -i + 1)
.

Let ρ ≡ ( 1 2 , . . . , 1 2 , 1 4 (1 -m)) and set ω m = m j=1 Γ (j/2) (2πi) m π m(m+1)/4 m! . ( 9
)
We shall use the notation

C m (ρ) = {s ∈ C m : Re(s) = -ρ} (10)
because this subset of C m arises frequently in the sequel, and we also define

d * s = ω m |c m (s)| -2 ds 1 • • • ds m (11)
since this measure is ubiquitous in our development.

Let M = {diag(±1, . . . , ±1)} be the collection of m × m diagonal matrices with entries ±1 on the diagonal; then M is a subgroup of K and is of order 2 m . By factorizing the Haar measure dk on K, it may be shown ( [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], p. 88) that there exists an invariant measure d k on the coset space K/M such that

k∈K/M d k = 1.
In stating the inversion formula for the Helgason-Fourier transform, we make particular use of the notation ( 9)- [START_REF] Muirhead | Aspects of Multivariate Statistics[END_REF]. The inversion formula then is the following result.

For f ∈ C ∞ c (P m ) and w ∈ P m ,

f (w) = C m (ρ) k∈K/M Hf (s, k) p s (k ′ wk) d k d * s. (12)
In particular, if f ∈ C ∞ c (P m /K) then

f (w) = C m (ρ) f (s) h s (w) d * s. (13)
We refer to Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], p. 87 ff. for a detailed treatment of this inversion formula and many references to the literature. For f ∈ L 1 (P m ) and h ∈ L 1 (P m /K) we define f * h, the convolution of f and h, by

(f * h)(w) = Pm f (z) h(z -1/2 wz -1/2 ) d * z,
w ∈ P m . Thus, if f and h are the density functions of independent random matrices X ∈ P m and ε ∈ P m /K, respectively, then f * h is the density function of the composition X • ε.

For f ∈ C ∞ c (P m ) and h ∈ C ∞ c (P m /K), the convolution property of the Helgason-Fourier transform is that

H(f * h)(s, k) = Hf (s, k) h(s), (14)
s ∈ C m , k ∈ K; see Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], Theorem 1, p. 88.

Deconvolution density estimation on P m

On the group G, the deconvolution problem arises from the statistical model

Y L = X ε ( 15
)
where X is a random unobservable, ε is an independent random error, and Y is the observed random measurement. We assume that f ε , the density of ε, is known and K-biinvariant and that the unknown densities f X and f Y of X and Y , respectively, are K-invariant. Under the equivalence K\G ≃ P m , we have ε → ε, X → X, and Y → Y together with the identification

f ε ( g) = f ε ( g ′ g), f X ( g) = f X ( g ′ g) and f Y ( g) = f Y ( g ′ g), g ∈ G.
Since X and ε are independent then X and ε also are independent, and ( 15) implies that

Y L = ε 1/2 Xε 1/2 ; ( 16
) hence f Y = f X * f ε .
Applying to [START_REF] Schwartzman | Cross-subject comparison of principal diffusion direction maps[END_REF] the convolution property ( 14) of the Helgason-Fourier transform, we obtain

Hf Y (s, k) = Hf X (s, k) f ε (s), (17)
s ∈ C m , k ∈ K.
Given a random sample Y 1 , . . . , Y n from Y , we estimate the density function f X as follows. We form H n f Y , the empirical Helgason-Fourier transform,

H n f Y (s, k) = 1 n n ℓ=1 p s (k ′ Y ℓ k). (18)
Substituting [START_REF] Zhu | Statistical analysis of diffusion tensors in diffusion-weighted magnetic image resonance data (with discussion)[END_REF] in [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], together with the assumption that f ε (s) = 0, s ∈ C m , we obtain

H n f X (s, k) = H n f Y (s, k) f ε (s) , s ∈ C m , k ∈ K.
In analogy with classical Euclidean deconvolution, we introduce a smoothing parameter T = T (n) where T (n) → ∞ as n → ∞, and then we apply the Helgason-Fourier inversion formula (12) using a spectral cut-off.

We introduce the notation

C m (ρ, T ) = {s ∈ C m (ρ) : λ s < T }
where C m (ρ) is defined in [START_REF] Matthews | Magnetic resonance imaging of multiple sclerosis: new insights linking pathology to clinical evolution[END_REF]. As an estimator of the population density f X , we define the density estimator f n X given by

f n X (w) = C m (ρ,T ) k∈K/M H n f Y (s, k) f ε (s) p s ( k′ w k) d k d * s,
w ∈ P m . The estimator f n X will serve as our nonparametric deconvolution estimator of the density f X .

Let {a n } and {b n } be sequences of real numbers. We write a n ≪ b n to mean a n ≤ Cb n for some constant C > 0, as n → ∞ (the Vinogradov notation). We use the notation a n = o(b n ) to mean a n /b n → 0, as n → ∞. We also write a n ≍ b n if both a n ≪ b n and b n ≪ a n ; and we write a n ∼ b n if a n /b n → 1 as n → ∞.

For technical reasons, we will also assume the moment condition

Pm |w 1 | -1 • • • |w m-1 | -1 |w| (m-1)/2 f Y (w) d * w < ∞, (19)
on the principal minors |y 1 |, . . . , |y m | of y ∈ P m . This assumption will be maintained throughout the rest of the paper. Estimation will proceed in the Sobolev class of functions,

H σ (P m ) = f ∈ C ∞ (P m ) : ∆ σ/2 f 2 < ∞ ,
where 2σ > dim P m = m(m + 1)/2, ∆ is the Laplacian on P m and, for f ∈ C ∞ (P m ),

f = Pm |f (w)| 2 d * w 1/2
denotes the L 2 (P m )-norm with respect to the invariant measure d * w. For Q > 0, we also define the bounded Sobolev class,

H σ (P m , Q) = f ∈ C ∞ (P m ) : ∆ σ/2 f 2 < Q , where 2σ > dim P m = m(m + 1)/2.
Theorem 1. Suppose there exists β ≥ 0 such that

| f ε (s)| -2 ≪ T β as T → ∞, for all s ∈ C m (ρ, T ). If f X ∈ H σ (P m , Q) and σ > 1 2 dim P m ≡ m(m + 1)/4 then, as n → ∞, E f n X -f X 2 ≪ n -2σ/(2σ+2β+dim Pm) . Corollary 1. Suppose the distribution of ε is concentrated at I m . If f X ∈ H σ (P m , Q) where σ > 1 2 dim P m then as n → ∞, E f n X -f X 2 ≪ n -2σ/(2σ+dim Pm) .
We shall also obtain a result for the situation in which the hypothesis in Theorem 1 is replaced by an exponential bound. In such a situation, we have the following result. Theorem 2. Suppose there exists β, γ > 0 such that

| f ε (s)| -2 ≪ exp(T β /γ), as T → ∞, for all s ∈ C m (ρ, T ). If f X ∈ H σ (P m , Q) with σ > 1 2 dim P m then, as n → ∞, E f n X -f X 2 ≪ (log n) -σ/β .
In this situation, we consider the special case in which

f ε (s) = exp(-γ -1 λ β s ), s ∈ C m ,
where γ > 0 is a scale parameter. Again by the inversion formula [START_REF] Pomara | White-matter abnormalities in HIV-1 infection: A diffusion tensor imaging study[END_REF], the underlying density function is

f ε (w) = C m (ρ) exp(-γ -1 λ β s ) h s (w) d * s,
w ∈ P m . The case in which β = 1 is particularly important and is called the heat or Gaussian kernel, since the latter is the fundamental solution to the heat equation on P m , see Terras [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], pp. 106-107. As a consequence, we obtain the following result.

Corollary 2. Suppose that f ε is Gaussian. If f X ∈ H σ (P m , Q) where σ > 1 2 dim P m then, as n → ∞,

E f n X -f X 2 ≪ (log n) -σ .
A case which is familiar in multivariate statistics is that of the Wishart distribution, W m (N, Σ), where Σ ∈ P m and N > m -1. For s = (s 1 , . . . , s m ) ∈ C m define the multivariate gamma function,

E f n X -f X 2 ≪ (log n) -2σ .

  For a, b ∈ C with Re(a), Re(b) > 0, let B(a, b) = Γ (a) Γ (b)/Γ (a + b) denote the well-known beta function, where Γ (•) is the classical gamma function. For s = (s 1 , . . . , s m ) ∈ C m , the Harish-Chandra c-function is

ΓTheorem 3 . 1 .

 31 m (s 1 , . . . , s m ) = π m(m-1)/4 m j=1 Γ s j + • • • + s m -1 2 (j -1) , where Re(s j + • • • + s m ) > (j -1)/2, j = 1, . . . , m. Relative to the invariant measure d * w in (6), the probability density function of the standard Wishart distribution W m (N, I m ) isf ε (w) = 1 Γ m (0, . . . , 0, N/2) m .We note that (20) is K-invariant and its Helgason-Fourier transform is well-known (Muirhead[START_REF] Muirhead | Aspects of Multivariate Statistics[END_REF], p. 248; Terras[START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications[END_REF], pp. 85-86),f ε (s) = Γ m ((0, . . . , 0, N/2) + s * ) Γ m (0, . . . , 0, N/2) h s ( 1 2 I m )wheres * = (s m-1 , s m-2 , . . . , s 2 , s 1 , -(s 1 + • • • + s m )). For N > m -1, the Wishart distribution W m (N, I m ) satisfies | f ε (s)| -2 ≪ exp πT 1/2 ,as T → ∞, where s ∈ C m (ρ, T ). Consequently we deduce the following result. Suppose that ε follows the Wishart distribution (20) with N > m -If f X ∈ H σ (P m , Q) with σ > dim P m /2 then, as n → ∞,

Research supported in part by the U.S. National Science Foundation, grants DMS-0112069 and DMS-0705210.


Vandana Mohan

Ganesh Sundaramoorthi

John Melonakos

Marc Niethammer

Marek Kubicki

Allen Tannenbaum

Tubular Surface Evolution for Segmentation of the Cingulum Bundle From DW-MRI

come   L'archive ouverte pluridisciplinaire

Introduction

In this work, we are interested in extracting a structure in the brain called the cingulum bundle (CB) from diffusion-weighted magnetic resonance imagery (DW-MRI) of the brain. DW-MRI is imagery that at each voxel indicates the diffusion of water molecules at each particular sampling direction in 3D space. Adding an extra dimension, directionality, to the data is necessary to discriminate our structure of interest -the cingulum bundle. The CB has recently become the subject of interest as an anatomical structure which may display quantifiable differences between schizophrenic and normal control populations, and studying it may aid in the diagnosis of schizophrenia [START_REF] Kubicki | Cingulate fasciculus integrity disruption in schizophrenia: a magnetic resonance diffusion tensor imaging study[END_REF][START_REF] Wang | Anterior Cingulum Abnormalities in Male Patients With Schizophrenia Determined Through Diffusion Tensor Imaging[END_REF].

The Cingulum Bundle is a thin, highly curved structure that consists of a collection of neural fibers, which are mostly disjoint possibly intersecting, roughly aligned and centered around a fiber. The collection of fibers approximately form a tube-like structure. The diffusion pattern in the CB varies in orientation and anisotropy smoothly along the structure, and it has a distinct diffusion pattern from surrounding areas of the brain (see Fig. 2 for a sagittal slice of the CB). The CB possesses a challenge to segment because of inhomogeneity of its diffusion pattern globally and the noisy nature of DW-MRI makes it difficult to detect edges separating the CB from the rest of the image.

There has been much research in detecting and characterizing neural connections between brain structures in DW-MRI. Early methods for detecting fibers, i.e., tractography, are based on streamlines where the fiber path is constructed by following the direction of the principal eigenvector of diffusion tensors from an initial seed point, e.g., [START_REF] Mori | Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF][START_REF] Conturo | Tracking neuronal fiber pathways in the living human brain[END_REF]. These methods have been shown to perform poorly in noisy situations and they often terminate prematurely before the fiber ending. To alleviate these problems, there has been a number of works, e.g., [START_REF] Perrin | Fiber tracking in q-ball fields using regularized particle trajectories[END_REF][START_REF] Parker | Estimating distributed anatomical brain connectivity using fast marching methods and diffusion tensor imaging[END_REF][START_REF] Friman | A bayesian approach for stochastic white matter tractography[END_REF][START_REF] Lenglet | A riemannian approach to diffusion tensor images segmentation[END_REF][START_REF]Control Theory and Fast Marching Techniques for Brain Connectivity Mapping[END_REF][START_REF] Zhang | Probabilistic fiber tracking using particle filtering[END_REF][START_REF] Cohen-Adad | Detection of multiple pathways in the spinal cord white matter using q-ball imaging[END_REF], where an optimal path, in some sense, is constructed from seed region(s). The procedure is repeated to detect all fibers of a bundle. These methods, however, do not explicitly provide a segmentation of the entire fiber bundle.

We are interested in segmenting the entire fiber bundle as a surface. Standard active surface techniques, e.g. [START_REF] Xu | Finding the brain cortex using fuzzy segmentation, isosurfaces, and deformable surface models[END_REF][START_REF] Yezzi | A geometric snake model for segmentation of medical imagery[END_REF], adapted to DW-MRI are typically difficult to segment the CB since the DW-MRI of the brain are extremely noisy and contain many local features that trap the active surface in unlikely configurations that are not representative of the CB. Indeed, the CB is difficult to segment without a shape prior favoring its thin tube-like geometry. Standard region-based techniques adapted to DW-MRI or DT-MRI, e.g. [START_REF] Lenglet | Dti segmentation by statistical surface evolution[END_REF], are generally not applicable to the segmentation of the CB since the statistics of the DW-MRI inside the CB cannot be described by a few global parameters (e.g. mean). The Mumford-Shah energy extended to DT-MRI, [START_REF] Wang | Dti segmentation using an information theoretic tensor dissimilarity measure[END_REF], which assumes piecewise smooth image data inside the surface, is applicable to the CB, but the technique needs a shape prior for the CB and is computationally costly since smooth functions must be determined at each update of the evolving surface. In [START_REF] Awate | Fuzzy nonparametric dti segmentation for robust cingulum-tract extraction[END_REF], the authors model the probability distribution of the CB and design an algorithm to classify voxels of the image, and the method could benefit greatly modeling the CB geometry. Noticing that standard region-based techniques are not applicable to the CB, an edge-based active surface method for segmenting the CB is considered by [START_REF] Melonakos | Finsler tractography for white matter connectivity analysis of the cingulum bundle[END_REF]. However, the method is sensitive to the noise in the DW-MRI, and the method does not incorporate the tube-like geometry into the segmentation. The work [START_REF] Melonakos | Locally-constrained region-based methods for dw-mri segmentation[END_REF] designs an energy on volumetric regions that incorporates "local region-based" information and a prior favoring regions that are close to an initially detected center-line curve. However, the energy is highly dependent on the correct placement of the detected center-line, which is often not exactly in the center of the CB. Moreover, the method does not enforce the tube-like geometry of the CB.

In this work, we explicitly model the CB as a tubular surface, and construct a general class of energies defined on these tubular surfaces. This enforces a tubular geometry during the segmentation process. Since the tubular surfaces we consider are determined by a center-line in 3D space and a radius function defined at each point of the center-line (see Fig. 1), the problem is reduced to optimizing an energy defined on 4D curves. This significantly reduces the computational cost of the optimization proce-dure when compared to an unconstrained surface optimization. Further, we show how to construct energies that are tailored to the varying nature of the diffusion pattern in the CB.

Our method is inspired by the work of [START_REF] Li | Vessels as 4-d curves: Global minimal 4-d paths to extract 3-d tubular surfaces and centerlines[END_REF] in which the authors model vessels as tubular regions formed by the union of spheres along a center-line. Energies are constructed on 4D curves that represent tubes, and these energies are globally minimized using the minimal path technique [START_REF] Cohen | Global minimum for active contour models: A minimal path approach[END_REF]. The energies we construct cannot be optimized using the minimal path technique since our energies are directionally dependent -they depend on the position of the 4D curve and its tangent. Moreover, for the energies we consider, we are not interested in a global minimum but rather certain local minima. As we shall see (in Section 4), the optimization of our energy of interest using gradient descent requires special consideration, and this interestingly ties to the metric structure on the space of 4D curves.

The Cingulum Bundle

The cingulum bundle is a 5-7 mm in diameter fiber bundle that interconnects all parts of the limbic system. It originates within the white matter of the temporal pole, and runs posterior and superior into the parietal lobe, then turns, forming a "ring-like belt" around the corpus callosum, into the frontal lobe, terminating anterior and inferior to the genu of the corpus callosum in the orbital-frontal cortex [START_REF] Schmahmann | Fiber Pathways of the Brain[END_REF]. Moreover, the CB consists of long, association fibers that directly connect temporal and frontal lobes, as well as shorter fibers radiating into their own gyri. The CB also includes most afferent and efferent cortical connections of cingulate cortex, including those of prefrontal, parietal and temporal areas, and the thalamostriatae bundle. In addition, lesion studies document a variety of neurobehavioral deficits resulting from a lesion located in this area, including akinetic mutism, apathy, transient motor aphasia, emotional disturbances, attentional deficits, motor activation, and memory deficits. Because of its involvement in executive control and emotional processing, the CB has been investigated in several clinical populations, including depression and schizophrenia. Previous studies, using DTI, in schizophrenia, demonstrated decrease of FA in anterior part of the cingulum bundle [START_REF] Kubicki | Cingulate fasciculus integrity disruption in schizophrenia: a magnetic resonance diffusion tensor imaging study[END_REF][START_REF] Wang | Anterior Cingulum Abnormalities in Male Patients With Schizophrenia Determined Through Diffusion Tensor Imaging[END_REF], at the same time pointing to the technical limitations restricting these investigations from following the entire fiber tract.

Proposed Framework

In this section, we model the cingulum bundle (CB) as a tubular surface in R 3 . We show that the tubular surface is completely determined by its center-line and the radius function of the discs along the center-line, and therefore, the tubular surface in R 3 can be effectively reduced to a curve in R 4 . We formulate a general class of energies directly on curves living in R 4 , and then observe that special consideration of the metric structure on curves is needed to optimize the energy.

Modeling the Cingulum Bundle (CB) as a Tubular Surface

We are interested in tubular surfaces since these surfaces naturally model the CB. These surfaces have the additional advantage that they may be represented as space curves thus significantly reducing the computational complexity of our algorithm.

S 1 × [0, 1] → R 3 (S 1 is [0, 2π]
with endpoints identified) as follows:

S(θ, u) = c(u) + r(u)[n 1 (u) cos θ + n 2 (u) sin θ] (1)
where n 1 , n 2 : [0, 1] → R 3 are normals to the curve c defined to be orthonormal, smooth, and such that the dot products c ′ (u) • n i (u) vanish. See Fig. 1 for an illustration of a tubular surface. The idea is simply that the tubular surface is represented as a collection of circles each of which lie in the plane perpendicular to the center line. Note that the surface in (1) may thus be identified with a 4D space curve, c : [0, 1] → R 4 , defined as a cross-product: c(u) = (c(u), r(u)) T .

(2)

Variational Approach for Detecting the Cingulum Bundle

We now define a general class of energy functionals defined directly on 4D curves (2) that when optimized result in the 4D curve that represents the CB from DW-MRI of the brain.

Let S 2 ⊂ R 3 denote the 2D sphere, which is to represent the set of all possible angular acquisition directions of the scanning device for DW-MRI. Let I : R 3 × S 2 → R + be the diffusion image. We are interested in weighted length functionals on 4D curves as energy functionals of interest. Indeed, let Ψ : R 4 × S 2 → R + (Ψ (x, r, v) ∈ R + ) be a weighting function, which we call the potential to be chosen, and define the energy as

E(c) = c Ψ (c(s), c ′ (s) |c ′ (s)| ) ds, c = (c, r) (3)
where

ds = |c ′ (u)| du = (r ′ (u)) 2 + |c ′ (u)| 2 du
is the arclength measure of the 4D curve, and c ′ (s)/|c ′ (s)| is the unit tangent to c, the center line. When (3) is minimized, the term ds penalizes the non-smoothness of the center line and the radius function.

The energy (3) is related to the length of a curve in a Finsler manifold [START_REF] Melonakos | Finsler active contours[END_REF].

The goal is to choose Ψ so that the energy is optimized by a c which determines a surface enclosing the diffusion pattern of the CB in the DW-MRI of the brain. The diffusion pattern in the cingulum varies in orientation and anisotropy across the length of the bundle, although locally similar (see Fig. 2), and that pattern differs from the pattern immediately outside the CB. This fact precludes the use of traditional region-based techniques adapted to DT-MRI since these techniques assume homogeneous statistics within the entire region enclosed by the surface, whereas we will assume homogeneity within local regions. In the next section, we show how Ψ may be chosen so that the energy can capture the varying diffusion pattern of the CB. The idea is to choose Ψ at a particular coordinate (x, r, p) to incorporate statistics of the DT-MRI local to the disc determined by (x, r, p) rather than using statistics global to the entire structure as in traditional region based methods.

Example Potentials, Ψ

In this section, we give two choices of Ψ that are meaningful for extracting the CB from DW-MRI, both based on local region-based statistics.

The first potential Ψ 1 at a coordinate (x, r, v) ∈ R 3 × R + × S 2 is constructed so as to be small when the mean diffusion profile inside the disc, D(x, r, v), differs greatly from the mean diffusion profile inside the annular region, D(x, αr, v)\D(x, r, v) where α > 1, outside D(x, r, p). This is given by the following expressions:

Ψ 1 (p, v) = 1 1 + µ D( p,v) -µ D((p,αr),v)\D( p,v) 2 (4)
where the µ's are means:

µ D( p,v) (v) = 1 r 2 D( p,v) I(x, v) dA(x) (5)
µ D((p,αr),v)\D( p,v) (v) = 1 (α 2 -1)r 2 D((p,αr),v)\D( p,v) I(x, v) dA(x), ( 6
)
where dA is the area element and • is a suitable norm on functions of the form f : S 2 → R + , e.g.,

f 1 -f 2 2 = S 2 |f 1 (v) -f 2 (v)| 2 dS(v), (7)
where dS is the surface area element. The energy corresponding to Ψ 1 is minimized. Another example potential is chosen such that the corresponding energy is related to a weighted surface area:

Ψ 2 (x, r, p) = r 2π 0 φ(x + rp ⊥ (θ)) dθ, and p ⊥ (θ) = n 1 cos θ + n 2 sin θ (8)
where n 1 , n 2 are orthonormal vectors perpendicular to p, and φ : R 3 → R + is large near the boundary of differing diffusion regions, e.g.,

φ(x) = 1 |B(x, R)| B(x,R) I(y, •) -µ B(x,R) (•) 2 dy (9)
where B(x, R) is the ball centered at x of chosen radius R, |B(x, R)| denotes the volume, and the norm is defined as in [START_REF] Friman | A bayesian approach for stochastic white matter tractography[END_REF]. For this choice of potential, we are interested in maximizing the corresponding energy. The objective is to initialize the tubular surface inside the CB, and then increase surface area until the surface reaches the boundary of differing diffusion patterns.

Energy Optimization

In this section, we construct a steepest descent flow to minimize the energy of interest [START_REF] Mori | Three dimensional tracking of axonal projections in the brain by magnetic resonance imaging[END_REF]. A steepest descent is considered since we are not necessarily interested in the global maximizer or minimizer; indeed, the energy corresponding to [START_REF] Lenglet | A riemannian approach to diffusion tensor images segmentation[END_REF] does not have a global maximizer. We begin with a tubular surface initialization (see Section 5 for the procedure), i.e., an initial 4D curve, and follow the gradient or its opposite depending on the whether we want to maximize or minimize the energy.

Gradient Descent: Fixed Endpoints

The standard technique for calculating the gradient of an energy defined on curves, which is based on a geometrized L 2 metric on the space of curves, cannot be applied to our energy of interest. This is because of the fact that when minimizing (3) using L 2 , Ψ must satisfy a certain positivity condition (see [START_REF] Melonakos | Finsler active contours[END_REF]) that we cannot guarantee for our choices of Ψ otherwise the gradient descent is ill-posed. Moreover, when maximizing (3), we are indeed maximizing a weighted length, which with respect to the standard geometrized L 2 curve metric, leads to an unstable reverse diffusion. As shown in [START_REF] Sundaramoorthi | New possibilities with Sobolev active contours[END_REF], such weighted length functionals may be optimized in a stable manner by moving in the gradient direction of the energy (3) with respect to a geometrized Sobolev metric:

Definition 1. Let c : [0, 1] → R 4 be such that c(0), c(1) are fixed. Let h, k : [0, 1] → R 4 be perturbations of c then h, k L 2 := 1 L c h(s) • k(s) ds, h, k Sob := L c h ′ (s) • k ′ (s)
ds, where L is the length of of the curve c, ds is arclength element of c, and the derivatives are with respect to the arclength parameter s.

It can be shown that the gradient of (3) with respect to the Sobolev metric above is

1 L ∇ Sob E(c) = K(Ψ p) + ∂ ŝK ( Ψv 1 + (r s/|c s|) 2 + Ψ cs ), (10)
where

K(f ) := L 0 K(•, s)f (s) ds, K(s 1 , s2 ) = 1 L s2 L (1 -s1 L ) 0 ≤ s2 ≤ s1 s1 L (1 -s2 L ) s1 ≤ s2 ≤ L . (11)
The expression has the additional numerical advantage that only first order derivatives are required in comparison to the standard L 2 gradient, which needs second order information, and simply cannot be used anyway since it results in an unstable flow. Note that as stated in [START_REF] Sundaramoorthi | New possibilities with Sobolev active contours[END_REF], the expression (10) may be computed efficiently in order N complexity, where N is the number of sample points of the curve.

Evolving Endpoints

In the previous subsection, we derived a gradient descent flow for (3) provided the endpoints of the 4D curve (i.e., the end cross sections of the tube) are fixed. We now describe how to evolve the endpoints so as to reduce the energy. This is useful for some choices of Ψ in (3), for example, Ψ 2 defined in [START_REF] Lenglet | A riemannian approach to diffusion tensor images segmentation[END_REF]. To determine the evolution of the endpoints, we compute the variation with respect to the endpoints. This results in

ct (0) = ∓ Ψv 1 + r s |c s| 2 ∓ Ψ cs , ct (1) = ± Ψv 1 + r s |c s| 2 ± Ψ cs , (12)
which will minimize/maximize the energy (depending on the sign chosen above). Therefore, the algorithm to reduce the energy is to alternatively evolve the endpoints by [START_REF] Xu | Finding the brain cortex using fuzzy segmentation, isosurfaces, and deformable surface models[END_REF] and then evolve the 4D curve by [START_REF] Zhang | Probabilistic fiber tracking using particle filtering[END_REF].

Experiments and Results

The algorithm was applied to DWI data of the brain from a data set that included schizophrenic and normal control subjects, with the DWI being acquired for 54 sampling directions. Results are included for the data from 2 subjects, and show the CB extracted for both the right and left bundles in each case.

In this paper, a perturbation of the anchor tract is used as the initial centerline curve and the smallest possible radius of 0.5 is used, with the surface essentially growing out from this initial radius. There are also other options to perform this initialization. Since we are given seed regions determined by an expert, an alternative initialization would be to connect the two seed regions with a streamline that passed above the Corpus Callosum (which is easy to segment). This initialization is being explored for future work.

The results included in this paper show the application of the proposed framework to the data sets, using the energy (3) using the potential Ψ 1 (4). Figure 2 shows slice-wise views of the CB segmentation results obtained from the proposed framework indicating the homogeneity of the discs within the captured volume. The Figures , 3, and4 show the tubular surface extracted by the proposed algorithm. It is to be noted that the surfaces are accurate while the boundaries shown in Figure 2 are the boundary locations rounded off to the grid points by the visualization process.

Conclusions and Future work

We have proposed a novel technique to extract the cingulum bundle, which is of interest in the medical community because of its ties to schizophrenia, from DW-MRI of the brain. Unlike other standard techniques for extracting fiber bundles in the brain, we are able to extract the entire bundle as a region at once rather than detecting individual fibers and then combining them to form the bundle, which is laborious and prone to errors, while also performing this as a curve evolution rather than a surface evolution thus avoiding the computational disadvantages of a levelset implementation. We have modeled the cingulum bundle as a tubular surface and constructed a variational approach to detect the optimum tubular surface from DW-MRI, which represents the CB. Tubular surfaces provide a natural and accurate shape prior for the cingulum bundle, and such a shape prior is necessary due to the noisy nature of the imagery and the fact that data is not very visible or highly corrupted in certain slices. As we have shown, the tubular surface can be represented as a 4D curve, and thus, we were able to significantly reduce the computational cost of the algorithm compared to extracting an arbitrary surface. The proposed model was shown to yield good segmentations of the Cingulum Bundle upon visual inspection; unfortunately, there is no expert ground truth data available since it is laborious to hand segment an entire volume and certain slices do not even display the CB diffusion pattern accurately.

In future work, the authors plan to explore different choices of Ψ in the energy functional(3), and explore smoothness terms for the tubes in the energy. We will also implement the evolution of endpoints, which will be useful since the initialization will have to only be a single seed point. Further, the use of the extracted Cingulum Bundles will be explored in population studies for the discrimination of Schizophrenia. We are also interested in applying the framework other tubular structures such as the Uncinate Fasciculus in the brain.

Fig. 1 .

 1 Fig. 1. Illustration of Tubular Surface model

Fig. 2 .

 2 Fig. 2. Selected slice-wise views of CB Segmentation results from proposed framework. The top row shows the DWI data and the bottom row shows the DWI data with the extracted surface rendered in 3D.

Fig. 3 .

 3 Fig. 3. CB Segmentation Results for Brain data set 1. Yellow shows the left CB and magenta shows the right CB.

Fig. 4 .

 4 Fig. 4. CB Segmentation Results for Brain data set 2. Yellow shows the left CB and magenta shows the right CB.

This work was supported in part by grants from NSF, AFOSR, ARO, MURI, as well as by a grant from NIH (NAC P41 RR-13218) through Brigham and Women's Hospital. This work is part of the National Alliance for Medical Image Computing (NAMIC), funded by the National Institutes of Health through the NIH Roadmap for Medical Research, Grant U54 EB005149. Information on the National Centers for Biomedical Computing can be obtained from http://nihroadmap.nih.gov/bioinformatics.


Thomas Albrecht
email: thomas.albrecht@unibas.ch

Reinhard Knothe
email: reinhard.knothe@unibas.ch

Thomas Vetter
email: thomas.vetter@unibas.ch

Modeling the Remaining Flexibility of Partially Fixed Statistical Shape Models

Statistical shape models are widely used to model the variability of biological shapes. They can be used to reconstruct missing information given partial or noisy data. In case of partial data, many different reconstructions are possible, and one is not only interested in a plausible reconstruction but also the remaining flexibility within the model and the reliability of the reconstruction. We present a method to model the remaining flexibility when some part of a statistical shape model is fixed. Using such a flexibility model, we can give answers to questions like: Does one half of a human femur bone determine the other half? or How much is the shape of a face determined by its contour?

Introduction

Statistical shape models are widely used in medical image analysis, computational anatomy, and computer vision to model the variability of biological shapes, see [START_REF] Fleute | Building a complete surface model from sparse data using statistical shape models: Application to computer assisted knee surgery[END_REF][START_REF] Rajamani | Statistical deformable bone models for robust 3d surface extrapolation from sparse data[END_REF][START_REF] Subsol | A scheme for automatically building threedimensional morphometric anatomical atlases: application to a skull atlas[END_REF][START_REF] Dam | Automatic shape model building based on principal geodesic analysis bootstrapping[END_REF][START_REF] Rao | Canonical Correlation Analysis of Sub-cortical Brain Structures Using Non-rigid Registration[END_REF][START_REF] Shen | An adaptive-focus statistical shape model for segmentation and shape modeling of 3-D brain structures[END_REF] for instance. The variability of a certain class of shapes is deduced from a representative set of example data from this class.

If the example data sets represent the class of shapes well, the model can be fitted to virtually any individual shape within the class. It is even possible to fit the model to partial data from an individual, [START_REF] Blanz | Reconstructing the complete 3d shape of faces from partial information[END_REF][START_REF] Keller | 3D Reconstruction of Human Faces from Occluding Contours[END_REF]. The missing remaining data is automatically reconstructed by the model, yielding the most plausible reconstruction.

However, there may be many other possible reconstructions which fit the partial data equally well. The partial data only determines a part of the model, while the rest of the model may remain flexible. This paper focuses on modeling this remaining flexibility of the partially determined model.

Keeping a part x b of the complete model x fixed, how much flexibility remains for the remaining part x a ? In principle, as PCA models can be statistically interpreted by a multivariate Gaussian distribution with probability density function p, we can model the remaining flexibility by the conditional distribution p(x a |x b ). We will see however, that there is no nontrivial conditional distribution if more components are kept fixed than there are degrees of freedom in the model, as is usually the case in models built from a small set of examples.

In this case, we propose a method which models the remaining flexibility of the variable points when the fixed points are allowed to move slightly instead of being completely fixed. The method leads to a generalized eigenvalue problem which can be solved efficiently.

The remaining paper is organized as follows. In Section 2 we will review the well-known concept of PCA-based statistical shape models. In Section 3 we will derive an expression of the conditional distribution p(x a |x b ) and investigate when this can actually be used. In Section 4, we will model the remaining flexibility when the fixed points are allowed to move slightly. In Section 5 we will use the proposed method to model the remaining flexibility of a shape model of the human femur bone when its distal part is fixed, and the remaining flexibility of a face model, when its contour is fixed.

Related Work. Virtually all cited papers deal with statistical shape models with applications in medicine or computer vision, [START_REF] Fleute | Building a complete surface model from sparse data using statistical shape models: Application to computer assisted knee surgery[END_REF][START_REF] Rajamani | Statistical deformable bone models for robust 3d surface extrapolation from sparse data[END_REF][START_REF] Subsol | A scheme for automatically building threedimensional morphometric anatomical atlases: application to a skull atlas[END_REF][START_REF] Dam | Automatic shape model building based on principal geodesic analysis bootstrapping[END_REF][START_REF] Rao | Canonical Correlation Analysis of Sub-cortical Brain Structures Using Non-rigid Registration[END_REF][START_REF] Shen | An adaptive-focus statistical shape model for segmentation and shape modeling of 3-D brain structures[END_REF][START_REF] Blanz | Reconstructing the complete 3d shape of faces from partial information[END_REF][START_REF] Keller | 3D Reconstruction of Human Faces from Occluding Contours[END_REF]. Most notably, [START_REF] Rao | Canonical Correlation Analysis of Sub-cortical Brain Structures Using Non-rigid Registration[END_REF] uses Canonical Correlation Analysis to model the correlation between different parts of a statistical model, computing model coefficients which maximize the correlation between these parts. However, they do not address our problem of modeling the remaining flexibility when one part of the model is fixed.

Principle Component Analysis

The 3D statistical shape models used in this paper are based on a Principle Component Analysis (PCA) of a set of training data comprised of n 3D surfaces. Each surface is represented by a triangular mesh with the same number m ∈ IN of vertices, which are stacked into a data vector x = (x 1 , y 1 , z 1 , ...x m , y m , z m ) T ∈ IR M with M = 3m. The training surfaces need to be in correspondence.

As in all PCA models, the samples are assumed to be i.i.d. samples drawn from a multivariate normal distribution that is approximated by the estimated multivariate normal distribution N (x, Σ). The mean is estimated by the arithmetic mean x = 1 n n i=1 x i of all samples. The covariance matrix Σ ∈ IR M ×M can be estimated from the mean-free data matrix X := [x 1x, . . . , x nx] ∈ IR M ×n as Σ = 1 n XX T . When X is decomposed with a (reduced) Singular Value Decomposition X = UWV T into the product of a column-orthonormal matrix U ∈ IR M ×n , a diagonal matrix W ∈ IR n×n , and an orthonormal matrix V ∈ IR n×n , the covariance matrix can be expressed as Σ = 1 n UW 2 U T . The columns u i of the matrix U are the eigenvectors of Σ. They are known as the principal components of the model and describe the main modes of variation of the training data. Their corresponding eigenvalues σ 2 i := 1 n w 2 i describe the variance of the model projected onto these eigenvalues. They are arranged according to size so that u 1 is the direction with maximal projected variance σ 2 1 . One individual x in the object class modeled by the PCA model can be identified by its coefficients α = (α 1 , . . . , α n ) T ∈ IR n :

x = x + n i=1 α i σ i u i = x + 1 √ n UWα. (1)
Under the assumption that the data x is distributed according to a multivariate normal distribution N (x, Σ), the coefficient vector α is distributed according to N (0, I n ).

In this paper two PCA models are used: A model of the femur bone, built from 21 CT scans of femur bones, which were hand-segmented and brought into correspondence with [START_REF] Dedner | Curvature guided level set registration using adaptive finite elements[END_REF], and a model of the human face built from 100 face scans that were brought in correspondence with a modification of [START_REF] Amberg | Optimal Step Nonrigid ICP Algorithms for Surface Registration[END_REF].

Conditional Distribution

In a PCA-based statistical model, we now wish to fix a certain number l ∈ IN of points in the model in order to investigate how flexible the model remains with respect to the remaining ml variable points. This means that we fix L = 3l components of the model vector x. Without loss of generality, we can assume that they are the last l components and x can be partitioned as x = (x a , x b ) T .

As x is distributed according to a multivariate normal distribution, the conditional distribution is a also a multivariate normal distribution N (µ a|b , Σ a|b ). Its mean and covariance can be calculated from x and Σ.

The matrix of principal components U ∈ IR M ×n can be partitioned according to x. For simplicity's sake, we define the matrix of the principal components scaled by the diagonal matrix W ∈ IR n×n as Q :

= 1 √ n UW ∈ IR M ×n . We have U = U a U b , Q = 1 √ n UW = 1 √ n U a W U b W = Q a Q b (2)
The covariance matrix Σ ∈ IR M ×M can be calculated from Q by Σ = 1 n UW 2 U T = QQ T . We can partition Σ as follows:

Σ = Σ aa Σ ab Σ ba Σ bb = Q a Q T a Q a Q T b Q b Q T a Q b Q T b . (3)
Note that it is only possible to calculate a nontrivial conditional distribution if the matrix Σ bb is invertible. A necessary condition for Σ bb to be invertible is that the number of fixed degrees of freedom L is less the number of training examples n of the statistical model. Let us, for the moment, assume that Σ bb is invertible and calculate the conditional distribution according. According to [START_REF] Bishop | Pattern recognition and machine learning[END_REF], the covariance matrix Σ a|b of the can be expressed as:

Σ a|b = Σ aa -Σ ab Σ -1 bb Σ ba (4) = Q a Q T a -Q a Q T b (Q b Q T b ) -1 Q b Q T a (5) = Q a I n -Q T b (Q b Q T b ) -1 Q b Q T a , (6)
where I n denotes the (n × n) identity matrix. While the original matrix U is column-orthonormal, the sub-matrix U b defined in Equation ( 2) in not. Therefore, we perform an additional a singular value decomposition

Q b = Ūb W b V T b
with a column-orthonormal Ūb . Thus, the expression can be expanded to:

Σ a|b = Q a I n -V b W b ŪT b Ūb W -2 b ŪT b Ūb W b V b Q T a (7) = Q a I n -V T b V b Q T a . (8)
Where, in the last step, we have used ŪT b Ūb = I n . Similarly, the mean µ a|b is given as:

µ a|b = xa + Σ ab Σ -1 bb (x b -xb ) (9) = xa + Q a V b W -1 b ŪT b (x b -xb ), (10)
with x = (x a , xb ) T .
Provided that Σ bb is invertible, x a can be reconstructed from x b as µ a|b . This is the reconstruction presented in [START_REF] Blanz | Reconstructing the complete 3d shape of faces from partial information[END_REF] as the maximum a posteriori reconstruction. All other reconstructions which fit x b are modeled by N (µ a|b , Σ a|b ). The more flexibility this distribution allows, the less reliable the reconstruction by µ a|b is.

For each of these reconstructions, the fixed values x b are matched equally well. Indeed, if we take a closer look at the expression for the covariance matrix Σ a|b in Equation ( 8) we notice that the inner part (

I n -V T b V b ) is a projection onto the orthogonal complement of the column space of V b . As V b is the "input matrix" of the SVD of Q b , this projection is a projection onto the kernel of Q b .
This means that the distribution N (µ a|b , Σ a|b ) models only linear combinations Q a α of the scaled principal components for which α is in the kernel of Q b , i.e. Q b α = 0. So for all α ∼ N (µ a|b , Σ a|b ), the deformation of the fixed points is zero.

Invertibility of Σ bb . So far, we have assumed that

Σ bb = Q b Q T b ∈ IR L×L is in- vertible, i.e. rank(Σ bb ) = L. However, if Q b has
less than L linearly independent columns we have rank(Σ bb ) < L. In particular, if the statistical model is built from less than L examples or less than L principal components are used, the rank of Σ bb will be less than L and the above calculations involving an inversion of Σ bb are not valid. In this case there is no nontrivial conditional distribution p(x a , x b ).

In this case the mean µ a|b can still be approximated by using the pseudoinverse of Σ bb instead of its inverse, [START_REF] Blanz | Reconstructing the complete 3d shape of faces from partial information[END_REF]. However, if we try to use the pseudoinverse for calculating the covariance matrix according to Equation (4), we get

Σ a|b = 0.
The problem is that the conditional distribution N (µ a|b , Σ a|b ) models only coefficient vectors from the kernel of Q b . Without a nontrivial kernel of Q b , it is not possible to calculate a nontrivial conditional distribution. Therefore, we that instead of looking for coefficient vectors α ∈ ker Q b , i.e. Q b α = 0, we will look for coefficients for which Q b α is small. This means that we relax the constraint of keeping the fixed points completely fixed to allowing them to move slightly.

It also has to be noted that even if there is a nontrivial kernel of Q b , it will most likely still be of interest to consider deformations which allow slight changes of the fixed points x b when considering the remaining flexibility of the model.

Generalized Eigenvalues

The aim is to model the flexibility of the variable points x a when the fixed points x b are allowed to move slightly. The deformations are given as Q a α and Q b α. A measure for the change caused by these deformations is the squared Euclidean norm of these vectors:

Q a α 2 = α T Q T a Q a α, b α 2 = α T Q T Q b α. ( 11
)
we are first and foremost interested in the coefficients α which change variable points x a as much as possible, we can formulate our aim as a constrained maximization problem:

max α∈IR k α T Q T a Q a α (12) subject to α T Q T b Q b α = c, (13)
where c ∈ IR + quantifies the amount of change allowed in the fixed coefficients. Introducing a Lagrangian multiplier λ and differentiating with respect to α leads to the generalized eigenvalue problem:

Q T a Q a α = λ Q T b Q b α. (14)
Both matrices Q T a Q a and Q T b Q b are positive definite and symmetric. The generalized eigenvalue problem can be solved efficiently with standard software (LAPACK, MATLAB), yielding a set of generalized eigenvectors {α 1 , . . . , α k } arranged according to the size of their corresponding generalized eigenvalues {λ 1 , . . . , λ k }. The eigenvectors are scaled so that

Q b α i 2 = α T i Q T b Q b α i = 1 for i = 1, . . . , k.
If we pre-multiply Equation ( 14) by α T , we see that for an eigenvector α i and its eigenvalue λ i , we have

Q a α i 2 = λ i Q b α i 2 . ( 15
)
This means that, measured in the squared Euclidean norm, the deformation determined by the coefficient vector α i ∈ IR k changes the variable model points λ i times as much the fixed ones. Therefore, the eigenvector α 1 corresponding to the largest eigenvalue λ 1 is the coefficient vector which causes the largest change on the variable points (with a squared Euclidean norm of λ 1 ), changing the fixed points only slightly (with a deformation with a squared Euclidean norm of 1). The last eigenvectors change the fixed points more than the variable ones. In

-2.2 √ l α1 +2.2 √ l α1 0 -2.2 √ l α2 +2.2 √ l α2
Fig. 1. In a statistical shape model of the human femur bone, the two first flexibility components α1 and α2 model those deformations that change the proximal part (gray) as much as possible, while changing the distal part (colored) as little as possible.

fact, exchanging fixed and points leads to the same eigenvalues, only inverted.

We call the coefficient vectors α i the flexibility components. α 1 is the first flexibility component. Just like the principal components of the PCA model, they can be scaled and added together to achieve different deformations from the mean. When a flexibility component α i is scaled by √ r it causes a deformation of squared Euclidean norm r of the fixed points and of λ i r of the variable points. The larger the eigenvalues λ i are, the more flexibility remains for the fixed points and the less reliable a reconstruction of x a from x b is considered. In order to meat the constraint from Equation (13), the vectors have to be scaled by √ c. The squared Euclidean norm is not a very intuitive measure of the deformation as it implicitly depends on the number l of fixed points. A more intuitive measure is the mean squared norm, averaged over all l fixed points:

1 l l k=1 (Q b α i ) 2 k = 1 l Q b α i 2 . ( 16
)
Therefore, if we scale the coefficient vectors α i by √ l, and the coordinates of the model are given in millimeters, the squared deformation of one model point is 1 millimeter on average. For the examples shown in Figures 4 and1, the deformation of the fixed vertices (blue) is compared with the flexibility of the full shape (red). The first flexibility components strongly affect the full shape, while they hardly alter the fixed vertices.

Regularization. The flexibility components α i are calculated as generalized eigenvectors. As outlined in Section 2, the coefficients α of the PCA model are distributed according to a multivariate normal distribution N (0, I n ). Without regularization, it is possible that some of the entries of α i can be extremely large. In terms of the distribution N (0, I n ), such extremely large values are very unlikely and in practice they cause unnatural deformations of the modeled shape. Up to a normalizing factor, the probability of a deformation caused by α can be calculated by e -1 2 α 2 . By minimizing α 2 = α T α, the probability of the deformation is maximized. Therefore, we propose replacing the constraint in Equation (13) by the following constraint:

α T Q T b Q b α + η α T α = c, (17)
a regularizing parameter η ∈ IR + . In this way, we limit not only the deformation of the fixed points given by Q b α 2 but also the improbability of the deformation, given by α 2 . The corresponding generalized eigenvalue problem is a regularized version of Equation ( 14) and is given by:

Q T a Q a α = λ (Q T b Q b + ηI n )α. (18)
This system no longer admits extremely large values in the generalized eigenvectors α i . The parameter η controls the balance between the original problem and the regularizing effect of allowing only probable shapes.

Note that when we deform not the mean but a certain individual x, which is defined by coefficients β, the deformation coefficients α are not distributed according to N (0, I n ) but according to N (-β, I n ). Therefore it would make more sense to minimize α + β 2 = (α + β) T (α + β). However, it is not obvious how to include such a constraint into a generalized eigenvalue problem like Equation ( 18). Therefore, we content ourselves with the proposed regularization, which penalizes large deformations, even though its statistical motivation is perfectible.

Nontrivial kernel of Q b . We have introduced the generalized eigenvalue problem Equation ( 14) in order to model the remaining flexibility of the variable points when it is not possible to calculate the conditional distribution p(x a |x b ), which is the case if and only if Q b has only the trivial kernel ker

Q b = {0}.
But what happens to the generalized eigenvalue problem if Q b does have a nontrivial kernel? In this case, we can split up the space of all coefficients α into the kernel of Q b and its orthogonal complement. In the kernel, we have deformations which do not change the fixed points at all, which can be modeled by the conditional distribution. The covariance matrix in the kernel is simply Q a Q T a , as in the kernel, the projection term from Equation ( 8), (I n -V T b V b ) = I n . In the complement, we can compute the generalized eigenvalue problem in order to additionally allow deformations which change the fixed points slightly.

From a practical point of view, when ker Q b = {0} and therefore Q T b Q b is singular, the aforementioned LAPACK or MATLAB routines return a basis for the kernel as generalized eigenvectors with eigenvalues infinity, as for these vectors we have Q T b Q b α i = 0. The remaining eigenvalues span the complement of the kernel and are computed as usual.

If we use the regularized form of the generalized eigenvalue problem Equation (18), there will be no infinite eigenvalues, as the matrix (Q

T b Q b + ηI n ) is always nonsingular. However, for vectors α ∈ ker Q b , we have (Q T b Q b + ηI n )α = ηα.
Therefore, on ker Q b , we effectively solve the eigenvalue problem:

α T Q T a Q a α = λ ηα, (19)
whose eigenvectors corresponding to the largest eigenvalues maximize Q a 2 . In this regularized case, there is no strict decomposition into ker Q b and its complement.

Experiments and Conclusion

We demonstrate the use of the flexibility models for two different scenarios with two different models, one for the shape of human faces and one for the human femur bone.

In the case of the femur bone, we are interested in determining how well the distal (bottom) part of the femur determines the proximal (top) part, in order to estimate how reliable the model can be for reconstruction of missing or injured parts of the bone. In Figure 1, the distal part of the bone is colored. It corresponds to the fixed points x b . From the generalized eigenvalue problem Equation (18), the flexibility components α i are calculated. The regularization parameter was chosen as η = 10. Figure 1 illustrates the effect of deforming the mean by ±2.2 √ l times the first two flexibility components α 1 and α 2 . Clearly, both the distal and the proximal part are changed, but the proximal part is much more heavily deformed. In Figure 2, the amount of deformation is plotted, measured in the mean and the maximum of the Euclidean norm of the deformation at each point. For the mean norm, the ratio is approximately equal to the corresponding eigenvalue. For the mean squared norm, which is not plotted here, it is of course exactly equal to this eigenvalue.

In case of the face model we are interested in the question: Given the vertices on the contour, how much is the shape of the face determined? Here, the fixed vertices x b are the vertices of the occluding contour in profile view, as shown in Figure (3 a)). This is the contour for one individual represented by its model coefficients β. Again, the flexibility components α i are computed using Equation (18), with η = 10. To visualize the result, the surface with the coefficients β ±1.1 √ lα i is shown in Figure 4. We see that the first two flexibility components heavily deform the model, while the vertices at the contour are almost fixed, see Figure 3 b,c). The resulting deformations measures are plotted in Figure 2. The ratio is much higher than in the femur case as less points are kept fixed and the model is built from many more examples, making it more expressive.

Conclusion.

We have introduced a way to model the remaining flexibility of a statistical shape model when a part of the model is kept as fixed as possible, even in the absence of a nontrivial conditional distribution. In future work, we will apply this technique in the fitting of statistical models to partial data.

  Fig. 2. For the examples shown in Figures4 and 1, the deformation of the fixed vertices (blue) is compared with the flexibility of the full shape (red). The first flexibility components strongly affect the full shape, while they hardly alter the fixed vertices.

Fig. 3 .

 3 Fig. 3. The vertices on the contour in profile view are fixed (a). The 1st (b) and 2nd (c) flexibility components (see Figure 4) don't change the contour in profile view.

Fig. 4 .

 4 Fig. 4. A statistical shape model of human faces is used to model the remaining flexibility with fixed contour. The vertices on the contour in profile view (see also Figure 3) are fixed. The flexibility components model the remaining variability within the model. The 1st/2nd flexibility components is shown in the 1st/2nd row.


N Kacem

S Baguet

R Dufour

S Hentz

Stability control of nonlinear micromechanical resonators under simultaneous primary and superharmonic resonances

Keywords: numbers: 85.85.+j, 05.45.-a, 46.40.Ff, 46.15.Ff

Fast effects of a slow excitation on the main resonance of a nonlinear micromechanical resonator are analytically and experimentally investigated. We show in particular how the bifurcation topology of an undesirable unstable behavior is modified when the resonator is simultaneously actuated at its primary and superharmonic resonances. A stabilization mechanism is proposed and demonstrated by increasing the superharmonic excitation.

Interests in nonlinear resonance have grown over the last decade leading to many applications, for example, in parametric amplification [3][4][5] and noise squeezing 2 eventually with bifurcation amplification 4 , switching of magnetization in single nanoparticles 5 , electromagnetic cancer detection through nonlinear resonance interaction [START_REF] Vedruccio | the Progress in Electromagnetics Research Symposium[END_REF] , spectroscopy to probe a threaded interface [START_REF] Riviére | [END_REF] , pull-in retarding under secondary resonances for MEMS resonators 8 as well as ultralow superharmonic resonance for functional nanowires 9 .

Here, unlike previous schemes under primary or single nonlinear resonance excitation, the mechanical structure is simultaneously driven under primary and superharmonic resonances. The effect of adding an extra harmonic excitation in a system has been investigated by very few authors, even theoretically. In the quasi-static case of a beam, it is shown that the axial buckling load can be increased by adding highfrequency excitation 10 . An analytical study of a Van Der Pol-Mathieu-Duffing resonator 11 showed a way to tune the nature of its nonlinear dynamic behavior by adding a fast harmonic excitation. Nevertheless, the fast effects of a slow harmonic excitation on the main resonance of a nonlinear system have neither been studied nor experimentally shown.

This paper analytically investigates a strategy of stability control based on an energy transfer within a single nonlinear system through simultaneous (primary and superharmonic) resonances. Furthermore, we employ this mechanism for the first time in micromechanical resonators to demonstrate experimentally its stability enhancement.

In order to understand simultaneous resonances, let us consider the following resonator equation, for a single degree-offreedom x:

ẍ + µ ẋ + ω n x + α 2 x 2 + α 3 x 3 +α 5 x 5 = ζ 1 cos (Ω) + ζ 2 cos (2Ω) (1)
The dot denotes derivation with respect to time, ω n the natural frequency, α 2 , α 3 , α 5 , the coefficients of the quadratic, cubic and quintic nonlinearities respectively, µ the constant damping coefficient, ζ 1 , ζ 2 the drive amplitudes of the first and the second harmonic of the drive frequency Ω.

The effect of each nonlinear term has been experimentally shown and explained, for both doubly clamped beams 12 and cantilevers 13 . The second-order term is due to the electrostatic force, but of course has no influence on the primary resonance like any even-order term. It is included here in the purpose of showing the effect of a superharmonic resonance. The third-order term is the well-known Duffing term, leading to bistable softening or hardening behavior depending on its sign. The fifth-order term is responsible for the amplitude response to become another multivalued (5 possible amplitudes) function of the frequency, a highly unstable so-called mixed behavior 14 . Micromechanical resonators are very commonly driven via the superposition of a static DC voltage V dc and a time varying AC voltage V ac . The electrostatic force being proportional to the square of the voltage, this scheme involves the presence of two harmonics Ω and 2Ω of Eq. ( 1), with amplitudes ζ 1 and ζ 2 proportional to V ac V dc and V 2 ac respectively. When Ω is tuned around ω n , the primary resonance at ω n is mainly actuated by the first harmonic Ω, the effect of the second harmonic being negligible. Experimentally, this corresponds to the "1f mode" typically used via a lock-in amplifier in MEMS and NEMS electrical characterizations. It has been previously shown that the onset of bistability of the Duffing behavior can be tuned via electrostatic negative nonlinearity: the third order term coefficient can then eventually be nulled so a stable and linear response is obtained 13,[START_REF] Nayfeh | Mook Nonlinear Oscillations[END_REF] . This is unfortunately true up to the onset of the mixed behavior only. Worse, investigation of the bifurcation topology in this "1f mode" showed the invariance of this onset, for a given resonator 14 : the onset of the mixed behavior sets a strict upper limit to the frequency stability of the resonator.

When Ω is tuned around ωn 2 , the resonant response at ω n is provoked by both the first and the second harmonic 2Ω. More precisely, in this "2f mode" scheme also commonly used, the 2Ω harmonic mainly generates the primary resonance at ω n , and the Ω harmonic actuates a superharmonic resonance at ω n . The response at ω n is thus made of simultaneous primary and superharmonic resonances. The latter one is generated via a "slow" excitation compared to the resonant frequency.

In order to investigate the fast effects of a slow excitation, we use the multiple scales method [START_REF] Nayfeh | Mook Nonlinear Oscillations[END_REF] to solve Eq. (1). A first-order uniform solution is sought in the form

x(t, ε) = x 0 (T 0 , T 1 ) + εx 1 (T 0 , T 1 ) + • • • (2)
where ε is the small nondimensional bookkeeping parameter, T 0 = t and T 1 = εt. Since both harmonics are linked by the drive frequency Ω, in order to analyze the nonlinear response under simultaneous primary and superharmonic resonances, we express the nearness of Ω to ωn 2 by introducing the detuning parameter σ according to

2Ω = ω n + εσ (3)
Substituting Eq. ( 2) into Eq. ( 1) and equating coefficients of like powers of ε yields for order ε 0

ζ 1 cos σT 1 + T 0 ω n 2 + ω 2 n x 0 + x 0 (2,0) = 0 (4)
for order ε 1

ζ 2 cos (σT 1 + T 0 ω n ) + x 2 0 α 2 + x 3 0 α 3 + x 5 0 α 5 +x 1 ω 2 n + cx 0 (1,0) + 2x 0 (1,1) + x 1 (2,0) = 0 (5)
where x

(j,k) i = ∂ k ∂T k 1 ∂ j x i ∂T j 0 .
The general solution of Eq. ( 4) can be written as

x 0 = X cos (ω n T 0 + Φ) - 4ζ 1 3ω 2 n cos ω n T 0 2 + σT 1 (6)
Equation ( 6) is then substituted into Eq. ( 5) and the trigonometric functions are expanded. The elimination of the secular terms yields two first order non-linear ordinary-differential equations which describe the amplitude X and phase β modulation of the response and allow for a stability analysis

Ẋ = - 4ǫζ 2 1 sin β 20Xα 5 ζ 2 1 cos β + 9α 2 ω 4 n 81ω 9 n - ǫµX 2 + O(ε 2 ) (7) β = - 2ǫζ 2 1 2α 2 cos β + 6α 3 X + 15α 5 X 3 9ω 5 n X +2σ - 40ǫα 5 ζ 4 1 {6 + cos(2β)} 81ω 9 n - 6ǫα 3 X 2 + 5ǫα 5 X 4 16ω n + O(ε 2 ) (8)
where β = 2σT 1 -Φ. The steady-state motions occur when Ẋ = β = 0, which corresponds to the singular points of Eqs. ( 7) and (8). Thus, the frequency-response equation under simultaneous primary and superharmonic excitations can be written and plotted parametrically with respect to the phase β.

The resonator design considered for analytical simulations is described in Fig. 1 Figure 2 shows the analytical nonlinear forced response function (FRF) of the beam resonator on its fundamental mode at ω n for Ω = ωn 2 (2f mode) for a fixed V ac = 0.3 V (hence a fixed 2Ω-excitation amplitude) and for different values of V dc . V dc = 0 means no superharmonic resonance and an increased value increases the superharmonic resonance effect. The latter exists when the second-order term becomes of significant value, and hence also when the higher-order terms do. Like Duffing primary resonance, bistability may be reached for simultaneous resonances at intermediate drive levels, and we will call its onset the critical amplitude. At these levels, the behavior is unchanged by simultaneous resonances: typical Duffing-like hardening (or softening) behaviors are observed, and the measured critical amplitude is equal to the primary resonance-only case: the slow excitation has little or no effect on the fast dynamics. When the drive level is still increased (and the fifth-order terms become non negligible), the question is wether or not simultaneous resonances response also exhibits the mixed behavior like primary-only resonance, and wether or not the onset of this mixed behavior can be tuned.

Figure 3 shows the response for different values of V dc (6, 8 and 10 V), and for a fixed V ac = 1 V , high enough to display a mixed behavior. The amplitude of the P bifurcation point, called the onset of the mixed behavior is shown. As opposed to the actuation under primary resonance only (see Fig. 4), the P point amplitude is shifted up by the increase of the AC voltage. This result opens up new ranges of resonator operation, in terms of control, stability and drive amplitude. Both slow (superharmonic) and fast (primary) dynamics are present in Eqs. ( 7) and ( 8), and modulate the amplitude X and the phase β of the system, due to terms proportional to cos β or sin β and to cos 2β respectively. Moreover, these terms are proportional either to α 5 or α 2 , coefficients of the quintic and quadratic nonlinearities producing the mixed behavior or internal resonance into the system respectively. The AC voltage sets the 2Ω excitation, while the DC voltage only amplifies the Ω excitation. The alteration of the bifurcation topology in the present simultaneous resonances configuration is thus only due to the increase of the superharmonic resonance.

The experimental investigation was carried out by placing the device in a vacuum chamber at ambient temperature. The discrepancy between the theoretical and measured resonance frequencies are explained by a residual stress of 15 M P a. A lock-in amplifier was used to measure the frequency response as shown in Fig. 1(c).

FIG. 4: Frequency responses showing mixed behaviors measured under primary resonance as well as under simultaneous primary and superharmonic excitations, the location of the different bifurcation points and the effect of the DC voltage on the P point.

For each 1f mode measurement, the traced voltage V out is the result of the vectorial substraction of the measurement with zero DC voltage from the AC+DC overall measurement, so the feedthrough voltage due to parasitics is gotten rid of and V out is only due to the mechanical motion of the resonator. For 2f mode measurements, the background due to feedthrough is already at a negligible level. Figure 4 shows six nonlinear mixed behavior peaks, dashed and solid curves corresponding to "1f mode" and "2f mode" respectively. The dimensionless voltage Ṽ = Vout V dc is plotted in Fig. 4. On the three 1f-mode curves where V ac = 0.5 V , the P point vertical location is unchanged with respect to the drive amplitude, proportional to V ac V dc . Incidentally, one can observe that the shape of the peak changes as the quality factor slightly varies with V dc due to ohmic losses [START_REF] Sazonova | A tunable carbon nanotube resonator[END_REF] . The onset of the mixed behavior is set only by the gap-to-width ratio of the device 14 .

The three 2f-mode curves were obtained for V ac = 1V , and a Vdc voltage increased from 6 to 10V. Since the Ωand 2Ω-excitation amplitudes are proportional to VacVdc and V 2 ac respectively, only the Ω-excitation amplitude was amplified. This increases significantly the quadratic nonlinearity, proportional to V 2 dc + 1 2 V 2 ac 12 , which in turn amplifies the superharmonic resonance. As a result, and due to the nonlinear interaction between the two harmonics, the slow dynamics transfers more energy into the resonance at ω n . The bifurcation topology is thus modified, and the onset of the mixed behavior in particular can be retarded as shown in Fig. 4: the slow Ωexcitation allows one to tune it; the simultaneous resonance shows the fast effect of a slow nonlinear resonance on the resonator bifurcation topology around the mixed behavior.

To conclude, generally the geometry of a resonator is designed to satisfy some goals in terms of performances (mass sensitivity, force sensitivity, frequency stability. . . ). These performances are limited by the signal to noise ratio, and hence the maximum vibration amplitude, usually set at the onset of bistability for stable operation and noise mixing rea-sons. One can use a capacitive gate to cancel out the third order nonlinear terms and operate linearly beyond the Duffing critical amplitude, which is possible in a certain range of DC and AC voltage, and width over gap ratio. This new stability domain is small and is jeopardized by the onset of the mixed behavior which cannot be tuned with only primary resonance. The use of simultaneous resonances (pri-mary+superharmonic) has been shown here to overcome this limitation, by stabilizing the dynamic behavior of the resonator when operated at high drives. In fact, the effect of the high order nonlinearities can be retarded by an energy transfer between fast and slow dynamics which enlarges the stability domain beyond the onset of the mixed behavior. As a result, a great gain in performance is expected via a large increase of the carrier power of the device. This should be particularly true with nanoscale devices which display reduced dynamic ranges or signal to noise ratio.

  (a). It uses 2 electrodes: electrode 1 for actuation and electrode 2 for sensing which allows for 2 port electrical measurements.All numerical simulations were carried out with the following set of parameters: l r = 200 µm, b = 4 µm, h = 2 µm, l a = 20 µm, g a = 1 µm, l d = 112 µm, g d = 750 nm, where h and b are the width and thickness of the microbeam, g a and g d are the actuation and sensing gaps. V ac and V dc were used for parametric investigations.

FIG. 1 :

 1 FIG. 1: (a): Optical microscope image of the device. 1 and 5: anchors, 2: microbeam, 3: actuation electrode, 4: sensing electrode. (b): SEM image of the device. (c): Connection layout for the electrical characterization.

FIG. 2 :

 2 FIG. 2: Analytical frequency responses under simultaneous primary and superharmonic excitations at intermediate drive levels Vac = 0.3V , showing hardening behavior, with varying V dc . Xmax is the displacement of the beam normalized by the gap g d at its middle point and is traced vs the detuning parameter σ to get rid of the frequency tuning due to change in V dc .