domain: math spi

text:
Giuseppe I Cirillo
email: gcirillo@ulg.ac.be

Alexandre Mauroy
email: a.mauroy@ulg.ac.be

Ludovic Renson
email: l.renson@bristol.ac.uk

Gaëtan Kerschen
email: g.kerschen@ulg.ac.be

Rodolphe Sepulchre
email: r.sepulchre@eng.cam.ac.uk

A spectral characterization of nonlinear normal modes

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

Introduction

Nonlinear normal modes (NNMs) of vibration can be regarded as the extension of the concept of linear normal modes to nonlinear systems. They have proved useful in a number of applications going from model reduction [START_REF] Touzé | Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures[END_REF] and the study of localization phenomena [START_REF] Vakakis | Normal modes and localization in nonlinear systems[END_REF] to system identification [START_REF] Peeters | Modal testing of nonlinear vibrating structures based on nonlinear normal modes: Experimental demonstration[END_REF][START_REF] Renson | Robust identification of backbone curves using control-based continuation[END_REF]. NNMs of undamped systems were first defined as a vibration of unison of the system [START_REF] Rosenberg | Normal modes of nonlinear dual-mode systems[END_REF]. The NNM definition was later generalized to encompass modal interactions [START_REF] Lee | Complicated dynamics of a linear oscillator with a light, essentially nonlinear attachment[END_REF]. Inspired by the center manifold technique, Shaw and Pierre extended the NNM concept to damped systems [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF]. They defined NNMs as two-dimensional invariant manifolds in phase space which are tangent to the planes formed by the linear modes at the equilibrium point. Such manifolds are invariant under the flow (i.e. trajectories initialized in the manifold do not leave the manifold), a property which extends the invariance of linear modes to nonlinear systems.

Operator-theoretic methods have a long and rich tradition in dynamical systems. In particular, a (nonlinear) system can be studied through the spectral properties of the associated (linear) Koopman operator. Koopman operators were first defined for Hamiltonian systems in [START_REF] Koopman | Hamiltonian systems and transformation in hilbert space[END_REF]. In the seminal work of Mezić [START_REF] Mezić | Spectral properties of dynamical systems, model reduction and decompositions[END_REF] the spectral properties of the operator were connected to geometrical features of ergodic (i.e. conservative) systems. More recently, this spectral approach was extended to the case of dissipative systems in the context of model reduction [START_REF] Mauroy | Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics[END_REF] and stability analysis [START_REF] Mauroy | Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control[END_REF]. The Koopman operator framework is also conducive to data analysis (see e.g. Dynamic Mode Decomposition [START_REF] Tu | On dynamic mode decomposition: Theory and applications[END_REF]) and was successfully used to capture coherent structures (called Koopman modes) or detect instabilities in various fields (e.g. fluid flows [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF], power grids [START_REF] Susuki | Nonlinear Koopman modes and power system stability assessment without models[END_REF], energy-efficient buildings [START_REF] Eisenhower | Decomposing building system data for model validation and analysis using the Koopman operator[END_REF]).

Our aim is to explore the connections that exist between NNMs and the Koopman operator framework. In this direction, a fundamental result reported by Mezić in [START_REF] Mezić | Spectral properties of dynamical systems, model reduction and decompositions[END_REF] is that normal modes of linear oscillations have their natural analogs, termed the Koopman modes, in the context of nonlinear dynamics. Specifically, it was shown that, for a hyperbolic fixed point, the dynamics can be linearized in the entire basin of attraction [START_REF] Lan | Linearization in the large of nonlinear systems and Koopman operator spectrum[END_REF] and that the coordinates for the linearization can be obtained through the eigenfunctions of the Koopman operator [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF]. This framework therefore seems very appropriate to extend the concept of linear normal modes to nonlinear systems. Indeed, we will demonstrate that there is a direct connection between NNMs and the spectral expansion of the Koopman operator. More precisely, the invariant manifolds are related to the zero level sets of some eigenfunctions of the operator, in direct analogy with the linear case.

Using this new characterization, we will show how issues related to the parametrization of the invariant manifold can be overcome. In the original work of Shaw and Pierre, a single pair of variables (either in physical or modal coordinates) was used to obtain a parametrization of the manifold, which can possibly hold only locally. For instance, the parametrization may fail when the invariant manifold presents a complex geometry with foldings, as shown in [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF][START_REF] Renson | An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems[END_REF]. In contrast, interpreting the concept of NNMs in the Koopman operator framework allows us to obtain a parametrization that is valid globally.

The paper is organized as follows. Sections 2 and 3 review the theory of NNMs and of the Koopman operator, respectively. In Section 4 the connection between the two frameworks is investigated, first in the linear case and then in the nonlinear case. We also derive a partial differential equation describing the geometry of the invariant manifold and an approximation of the manifold using Taylor series. The theoretical developments are illustrated with a twodegree-of-freedom mechanical system with cubic nonlinearity in Section 5. The conclusions of the present study are summarized in Section 6.

Introduction to NNMs

Consider the first-order equations of motion of an N-degree-of-freedom nonlinear me-

chanical system ẋi = y i , ẏi = f (x, y) i = 1, ..., N (1)
where x i represents a generalized coordinate (displacement or rotation), y i is the corresponding velocity, x = [x 1 , ..., x N ] and y = [y 1 , ..., y N ]. From this point on, we assume that system (1) is damped, (x, y) = (0, 0) is an equilibrium point, and the system linearized around (0, 0) has n/2 pairs of complex-conjugate eigenvalues with non-zero real parts (hyperbolic equilibrium point).

In the 1990s, Shaw and Pierre defined a NNM as a two-dimensional manifold in phase space [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF] that is invariant under the flow, i.e. trajectories with initial conditions in the manifold remain in it for all time. As such, the invariance property of linear normal modes was extended to nonlinear systems: a motion that starts on a mode stays on this mode, while the other modes remain quiescent. Inspired by the center manifold technique [START_REF] Carr | Applications of Centre Manifold Theory[END_REF], Shaw and Pierre mathematically described the two-dimensional invariant manifold using a pair of master coordinates, the other state-space coordinates being functionally related to the master pair. Choosing arbitrarily (x i , y i ) as master coordinates, the remaining coordinates follow the constraint equations

x j = X j (x i , y i ), y j = Y j (x i , y i ), j = 1, ..., N ; j = i. ( 2
)
In several other works (see, for instance, [START_REF] Pesheck | A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds[END_REF]), a linear change of coordinates was used to express Eq. ( 2) in modal space and describe the invariant manifold as a function of a linear modal coordinate and its corresponding modal velocity. A set of equations for the functions X j and Y j can be derived by substituting the time derivative of Eq. ( 2) into the equations of motion [START_REF] Touzé | Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures[END_REF]. This leads to a set of 2N -2 coupled partial differential equations (PDEs)

Y j (x i , y i ) = ∂X j (x i , y i ) ∂x i y i + ∂X j (x i , y i ) ∂y i f i , f j = ∂Y j (x i , y i ) ∂x i y i + ∂Y j (x i , y i ) ∂y i f i , (3)
with

f j = f j (x i , X, y i , Y), X = {X j (x i , y i )}, Y = {Y j (x i , y i )}
, and j = 1, ..., N ; j = i. The PDEs (3) do not depend on time and can be solved for X and Y. Finally, after a substitution of the solutions in the ordinary differential equations governing the master coordinates, the dynamics on the NNM is reduced to a single-degree-of-freedom (SDOF) nonlinear oscillator

dynamics ẋi = y i , ẏi = f i (x i , X, y i , Y). (4)
The PDEs (3) admit several solutions that correspond to the extension of the underlying linear normal modes. The constraint relations (X, Y) describe the geometry of a single NNM in phase space. For a linear system, the mode corresponds to a plane whose dependence in x i and y i is function of the inertia, stiffness and damping properties of the system. For a nonlinear system, the mode is no longer flat but curved due to nonlinear distortions. At the equilibrium point, the NNM is tangent to the plane formed by the linear mode of the underlying linear system. Figure 1 illustrates the in-phase NNM of a two-degree-of-freedom system with cubic nonlinearity. The surface shown in the figure corresponds to the position of the second DOF x 2 as a function of the master coordinates (x 1 , y 1 ). The solutions of (3) can be approximated in the form of a polynomial series expansion in x i and y i [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF], whose coefficients are the solutions of a set of algebraic equations. More recently, several numerical methods were proposed in the literature and applied to various conservative [START_REF] Pesheck | A new Galerkin-based approach for accurate non-linear normal modes through invariant manifolds[END_REF][START_REF] Blanc | On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems[END_REF] and nonconservative [START_REF] Renson | An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems[END_REF] systems, see [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF] for a review.

Figure 1: In-phase NNM of a two-degree-of-freedom system with cubic nonlinearity. The manifold was computed using the Taylor series expansion presented in Section 4.4 .

Either derived using state space or modal coordinates, Eq. (2-4) are based on the fundamental underlying assumption that the invariant manifold can be described as a function of the chosen pair of variables. However, the invariant manifold may present a complex geometry with several foldings that appear for a number of reasons including internal resonances, localization, and multiple fixed points [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF][START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF][START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF]. Such a manifold cannot generally be described using a simple parametrization that is linear in the state space coordinates. For instance, the invariant manifold shown in Figure 1 exhibits foldings, which prevent a parametrization in terms of the master coordinates (x 1 , y 1 ). In this case, a proper parametrization cannot be obtained as a solution of the PDEs (3).

In Ref. [START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF], the parametrization issue was addressed using a larger set of master coordinates. When considering n m pairs, the slave coordinates can be described as

x j = X j (x m , y m ), y j = Y j (x m , y m ), (5)
where (x m , y m ) are the vectors of master coordinates. After solving the corresponding manifold-governing PDEs, the dynamics is reduced to m coupled (nonlinear) oscillators. Although this method is elegant and proved effective on a nonlinear beam [START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF][START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF], it does not completely solve the intrinsic parametrization issue and may therefore fail in some regions of the phase space.

In this paper the parametrization issue is addressed using the relation that exists between NNMs and Koopman operator. Through a nonlinear change of coordinates, the dynamics on the NNM will be expressed as a single-degree-of-freedom linear oscillator. This novel approach bears strong resemblance with the normal form method [START_REF] Touzé | Nonlinear normal modes for damped geometrically nonlinear systems: Application to reduced-order modelling of harmonically forced structures[END_REF][START_REF] Neild | Nonlinear modal decomposition using normal form transformations[END_REF], whose main idea is to simplify the vector field using a series of coordinate changes. Unlike the normal form method which has only local validity, the Koopman operator framework provides a change of coordinates that is valid in the entire basin of attraction of a hyperbolic fixed point [START_REF] Lan | Linearization in the large of nonlinear systems and Koopman operator spectrum[END_REF]. Therefore, it is possible to think of methods based on Koopman operator as a generalization of methods based on normal form, with the notable difference that the system of coordinates obtained in this case is valid in the whole basin of attraction of a fixed point. By doing so the parametrization issue is overcome, and coupled PDEs describing the NNM in the entire basin of attraction of the equilibrium point can be obtained.

Koopman Operator

Motivation and definition

k a c k nl m 1 x 1 k b c m 2 x 2 k a c
Figure 2: 2DOF system with cubic nonlinearity.

The system depicted in Figure 2 is considered in order to illustrate the Koopman operator definition and properties. The system dynamics can be expressed through a system of ordinary differential equations (ODEs)

ẋ = f (x), x = (x 1 , x 2 , x 3 , x 4 ) T ; (6)
where x 1 and x 2 are, respectively, the positions of the first and second mass, and x 3 and x 4 the corresponding velocities. The general solution, or flow, of ( 6) is a function, φ φ φ(x 0 , t), that assigns a point in state space to an initial condition x 0 and a time t.

It is often interesting to analyze the dynamics of a system in terms of specific quantities such as, for instance, the relative displacement between the first and second mass g

1 (x) = x 2 -x 1 , the potential energy g 2 (x) = 1 2 (k a (x 2 1 + x 2 2 ) + k b (x 1 -x 2 ) 2 ) + 1 4 k nl x 4 1 or the kinetic energy g 3 (x) = 1 2 (m 1 x 2 3 + m 2 x 2 4
). All these quantities are functions from the state space to R (or more generally C) and are called observables. 5

Since state-space variables change with time, so do these observables. Given initial conditions x 0 , the time evolution of an observable is obtained through its composition with the flow. For instance, the potential energy varies according to V (x 0 , t) = g 2 (φ φ φ(x 0 , t)), while the kinetic energy is given by K(x 0 , t) = g 3 (φ φ φ(x 0 , t)). At every time instant t, K(x 0 , t) and V (x 0 , t) can still be considered as observables (i.e. functions from the state space to C) such that, generally speaking, the time evolution of an observable results in a family of observables given by its composition with the flow.

The Koopman operator theory concentrates on how arbitrary observables change with time, instead of how trajectories are organized in state space. Given a dynamical system

ẋ = f (x), x ∈ R n , (7)
and a vector space of observables F, the Koopman operator (or more rigorously, semigroup of operators) U t is defined as the composition of the flow with g(x) ∈ F:

U t g(x) = g(φ(x, t)). ( 8
)
Equation ( 8) considers the case of a scalar observable g. The definition can be extended to a vector-valued function g(x) by applying U t to each component g i (x) of the vector, provided that each vector component g i (x) is in F.

The fundamental property of Koopman operator is linearity. Recalling the example in Figure 2, the evolution of the total energy of the system, g 4 (x) = g 2 (x) + g 3 (x), is given by the sum of the evolutions of g 2 and g 3 , i.e. E(x 0 , t) = K(x 0 , t) + V (x 0 , t). More generally, applying U t to the linear combination of two observables g α (x) and g β (x) that belong to F yields

g γ (x) = a g α (x) + b g β (x) ( 9
)
U t g γ (x) = g γ (φ φ φ(x, t)) = a g α (φ φ φ(x, t)) + b g β (φ φ φ(x, t)) = = a U t g α (x) + b U t g β (x). ( 10
)
The linearity property of the Koopman operator is general and does not depend on the dynamical system. This property is very attractive, especially in the context of vibration analysis, where eigenvalues and eigenvectors form the basis of vibration modes. The Koopman operator allows such tools to be used with nonlinear systems, but the price to pay is that infinite dimensional vector spaces have to be considered.

Spectral properties

From this point on the space of observables F is fixed. Although the only requirement is invariance with respect to the action of U t (given g ∈ F, U t g ∈ F for every t), different choices of F result in different spectral properties. In this work, F is taken as the space of C 1 functions that are analytic in a neighborhood of the origin. Moreover, only systems with a stable hyperbolic fixed point (assumed to be the origin) are considered, and our investigations are restricted to the basin of attraction of that point.

Based on the above assumptions, this section presents the Koopman operator spectral properties that are essential to this paper. A more comprehensive introduction can be found in Ref. [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF].

Eigenfunctions

Since U t is linear, it is possible to consider its eigenfunctions: those observables s(x) ∈ F whose evolution in time is given by

U t s(x) = s(φ φ φ(t, x)) = exp(λt)s(x) ∀t, (11)
in which λ ∈ C is the corresponding eigenvalue.

As an example consider a simple one-dimensional system

ẋ = -ηx, η > 0, (12)
with corresponding flow φ(x, t) = x exp(-ηt). Applying the Koopman operator to the observable s 1 (x) = x yields

U t s 1 (x) = s 1 (φ(x, t)) = x exp(-ηt) = exp(-ηt)s 1 (x), (13)
thus, the identity function s 1 (x) = x is an eigenfunction with eigenvalue -η. Now consider a general linear system defined by a matrix A. The eigenvalues of A are also part of the spectrum of the Koopman operator. Moreover, this property extends to nonlinear systems: if A is the Jacobian matrix evaluated at the origin of the vector field defined in [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF], that is

A = ∂f ∂x x=0 , (14)
then its eigenvalues λ k are part of the spectrum of the Koopman operator [START_REF] Mauroy | Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control[END_REF]. In fact, if A admits n linearly independent eigenvectors, its eigenvalues are associated with n independent eigenfunctions s k (x) that are defined in the entire basin of attraction of the equilibrium. However, the spectrum of U t is much richer than the spectrum of A: it always contains an infinity of eigenvalues. For instance, in the case of system ( 12), the observable s 2 (x) = x 2 is also an eigenfunction:

U t s 2 (x) = (φ(x, t)) 2 = (x exp(-ηt)) 2 = x 2 exp(-2ηt) = exp(-2ηt)s 2 (x). (15)
A similar reasoning can be applied to each integer power n of x, so that -nη and s n (x) = s n 1 (x) = x n are eigenvalues and eigenfunctions of system [START_REF] Tu | On dynamic mode decomposition: Theory and applications[END_REF], respectively. It turns out that this construction can be applied to general systems. Let s 1 (x) and s 2 (x) be two eigenfunctions:

U t s 1 (x) = exp(λ 1 t)s 1 (x) U t s 2 (x) = exp(λ 2 t)s 2 (x). ( 16
)
Applying the Koopman operator to their product s 1,2 (x) = s 1 (x)s 2 (x) yields

U t s 1,2 (x) = s 1,2 (φ φ φ(x, t)) = s 1 (φ φ φ(x, t))s 2 (φ φ φ(x, t)) = = exp(λ 1 t)s 1 (x) exp(λ 2 t)s 2 (x) = exp((λ 1 + λ 2 )t)s 1 (x)s 2 (x). ( 17
)
Thus, the sum of two eigenvalues is also an eigenvalue, and equivalently the product of two eigenfunctions is an eigenfunction. Using [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF] it is possible to construct an entire system of eigenfunctions. Considering the basis eigenfunctions s k associated with the eigenvalues λ k of the Jacobian matrix, we can obtain an infinity of eigenfuctions and eigenvalues:

s k 1 ,...,kn (x) = s k 1 1 (x) . . . s kn n (x), λ k 1 ,...,kn = k 1 λ 1 + . . . + k n λ n . ( 18
)
In the following a particular property of the eigenfunctions will be central to establish a connection to NNMs defined in Section 2. Let s(x) be an eigenfunction and λ the corresponding eigenvalue; then, its zero level set {x ∈ R n | s(x) = 0} is invariant for the dynamics (7) [START_REF] Mauroy | Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control[END_REF]. Indeed, given an initial condition x 0 satisfying

s(x 0 ) = 0, ( 19
)
the corresponding trajectory satisfies

s(φ φ φ(t, x 0 )) = exp(λt)s(x 0 ) = 0. ( 20
)
If the eigenvalue associated with s(x) is real, then so is s(x), and (19) corresponds to one equality. If λ has a nonzero imaginary part, s(x) cannot be purely real, so that ( 19) is a complex equation yielding two equalities Re(s(x)) = 0 and Im(s(x)) = 0. In the latter case, the complex conjugate λ is also an eigenvalue with eigenfunction s k (x), so that [START_REF] Renson | An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems[END_REF] corresponds to the zero level set of two (complex conjugate) eigenfunctions.

Spectral expansion

One of the main tools in the analysis of linear systems is the expansion of the dynamics in terms of linear modes. For a general linear system

ẋ = Ax, x ∈ R n , ( 21
)
with A diagonalizable, the flow can be written as

φ φ φ(x, t) = n k=1 s k (x)v k exp(λ k t), ( 22
)
where v k is a (possibly complex) eigenvector of A, λ k is the corresponding eigenvalue, and s k (x) is the projection of x onto the eigenspace generated by v k . Note the use of the same notation for the projections of x and the Koopman operator eigenfunctions; these two concepts will be linked in Section 4.1.

The Koopman operator theory allows the decomposition [START_REF] Blanc | On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems[END_REF] to be generalized to nonlinear systems [START_REF] Mezić | Analysis of fluid flows via spectral properties of the Koopman operator[END_REF]. Let g be an observable that belongs to the span of the eigenfunctions [START_REF] Jiang | The construction of non-linear normal modes for systems with internal resonance[END_REF]. The observable can be expressed as a (possibly infinite) sum of eigenfunctions

g(x) = (k 1 ,...,kn)∈N n k 1 +...+kn>0 C k 1 ,...,kn s k 1 1 (x) . . . s kn n (x). ( 23
)
Applying U t to [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF], the time evolution of the observable is given by

U t g(x) = (k 1 ,...,kn)∈N n k 1 +...+kn>0 C k 1 ,...,kn s k 1 1 (x) . . . s kn n (x) exp((k 1 λ 1 + • • • + k n λ n )t). ( 24
)
The coefficients C k 1 ,...,kn are the so-called Koopman modes. Each of them is the projection of the observable g onto the eigenspace generated by the corresponding eigenfunction s k 1 1 . . . s kn n . If a vector-valued observable g : R n → C d is considered, the Koopman modes in [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF][START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF] become vector-valued modes C k 1 ,...,kn and the scalar eigenfunctions remain unchanged.

As an example, consider for system (12) the observable g 1 (x) = x + x 2 . In the previous section it was shown that the functions s k 1 (x) = x k are eigenfunctions. Since g 1 (x) is a finite sum of eigenfunctions, we directly find the corresponding Koopman modes:

C 1 = 1, C 2 = 1, C k = 0, k ≥ 3.
A less trivial example of spectral decomposition is given by the function g 2 (x) = exp(x) -1. In this case, the Taylor series

g 2 (x) = exp(x) -1 = ∞ k=1 1 k! x k = ∞ k=1 1 k! s k 1 (x), (25)
shows that the Koopman modes are C k = 1 k! . In both cases [START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF] can be used to obtain the time evolution of the observables.

It is clear that [START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF] bears strong resemblance with the linear mode expansion [START_REF] Blanc | On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems[END_REF]. To obtain the same decomposition of the flow as in [START_REF] Blanc | On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems[END_REF], it is necessary to find an observable whose evolution corresponds to the flow. This is achieved using the identity function id(x) = x. In a neighborhood of the origin, id(x) can be decomposed as sum of eigenfunctions [START_REF] Mauroy | Isostables, isochrons, and Koopman spectrum for the action-angle representation of stable fixed point dynamics[END_REF]

id(x) = x = (k 1 ,...,kn)∈N n k 1 +...+kn>0 v k 1 ,...,kn s k 1 1 (x) . . . s kn n (x). (26)
Using [START_REF] Neild | Nonlinear modal decomposition using normal form transformations[END_REF] together with [START_REF] Boivin | Non-linear modal analysis of structural systems featuring internal resonances[END_REF], a decomposition of the flow in terms of eigenfunctions is obtained

U t id(x) = id(φ φ φ(x, t)) = φ φ φ(t, x) = (k 1 ,...,kn)∈N n k 1 +...+kn>0 v k 1 ,...,kn s k 1 1 (x) . . . s kn n (x) exp((k 1 λ 1 + . . . + k n λ n )t). ( 27
)
The decomposition [START_REF] Arnol'd | Geometrical methods in the theory of ordinary differential equations[END_REF] generalizes [START_REF] Blanc | On the numerical computation of nonlinear normal modes for reduced-order modelling of conservative vibratory systems[END_REF] to nonlinear systems. This infinite series can be truncated to a finite number of terms to approximate the system trajectories. The Koopman modes v k 1 ,...,kn associated with the identity function can be computed by solving a sequence of linear systems; the details of the calculation are reported in appendix A. For linear systems, Eq. ( 27) coincides with the linear normal mode decomposition, as shown in Section 4.1.

Koopman Operator and Vibration Modes

The linear case

Consider the linear system ẋ = Ax, x ∈ R n , [START_REF] Haller | Nonlinear normal modes and spectral submanifolds: Existence, uniqueness and use in model reduction[END_REF] with A diagonalizable, and the corresponding expansion in linear modes given by

φ φ φ(x, t) = n k=1 s k (x)v k exp(λ k t). ( 29
)
The terms s k (x) are the projections onto the corresponding eigenspace; they can be written in the form

s k (x) = w * k x, (30)
where * denotes the conjugate transpose and w k are the eigenvectors of A * , normalized so that

w * k v h = 1, k = h 0, k = h . ( 31
)
Applying the Koopman operator to s k (x), using the previous property, yields

U t s k (x) = s k (φ φ φ(x, t)) = w * k n i=1 s i (x)v i exp(λ i t) = s k (x) exp(λ k t). ( 32
)
Thus, the projections s k (x) are eigenfunctions of U t (which explains the notation used for them). It follows that the expansions of the flow through linear modes and Koopman operator eigenfunctions are identical in the linear case. Also, the general decomposition [START_REF] Arnol'd | Geometrical methods in the theory of ordinary differential equations[END_REF] appears as the natural extension of linear normal mode decomposition to nonlinear systems. The relation between vibration modes and spectral expansion is well-known in the linear case: given a pair of complex conjugate eigenvalues λ 1 and λ 2 , the motion on the linear mode is obtained when (29) reduces to

s 1 (x)v 1 exp(λ 1 t) + s 2 (x)v 2 exp(λ 2 t), (33)
i.e. when the following conditions on the eigenfunctions are verified:

s k (x) = 0, k = 3, . . . , n. ( 34
)
Thus, the plane defining the linear mode is determined by Eq. (34). Given a pair of complex conjugate eigenvalues, the corresponding linear mode is the intersection of the zero level set of all the other eigenfunctions.

In the next two sections (34) is generalized to nonlinear systems, and it is shown that this generalization is equivalent to define NNMs in terms of invariant manifolds.

Nonlinear coordinate transformation

The properties of Koopman operator in the general case of a nonlinear system with hyperbolic fixed point can be studied through linearization techniques [START_REF] Lan | Linearization in the large of nonlinear systems and Koopman operator spectrum[END_REF][START_REF] Mauroy | Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control[END_REF]. The results of [START_REF] Lan | Linearization in the large of nonlinear systems and Koopman operator spectrum[END_REF] are particularly useful, as they allow a local smooth change of coordinates (or diffeomorphism) h(x) to be extended to the entire basin of attraction of the fixed point. Assuming that the eigenvalues are nonresonant (i.e. relations of the type i n i λ i = 0 with n i ∈ N and i n i > 0 are not possible), there exists in a neighborhood of the origin an analytic diffeomorphism of the form [START_REF] Arnol'd | Geometrical methods in the theory of ordinary differential equations[END_REF]

z = x + Z(x), Z(0 0 0) = 0 0 0, ∂Z ∂x x=0 0 0 = 0 0 0, ( 35
) such that ż = Az. ( 36
)
Note that normal form methods aim at approximating the inverse transformation x = z + X(z) through its Taylor series.

According to the results and methods given in [START_REF] Lan | Linearization in the large of nonlinear systems and Koopman operator spectrum[END_REF], there exists a linearizing global C 1 diffeomorphism h(x), which is analytic when restricted to a neighborhood of the origin and characterizes the eigenfunctions of the Koopman operator. In particular, the eigenfunction corresponding to the eigenvalue λ k of A is given by

s k (x) = w * k h(x), (37)
where w k is a left eigenvector of A (here and in the following A is assumed diagonalizable). Moreover, it follows from (35) that

∂h ∂x x=0 0 0 = I, (38)
so that

∂s k ∂x x=0 0 0 = w * k , (39)
i.e. the gradient of the eigenfunction at the origin is a left eigenvector of A. Note that the eigenfunctions s k (x) = w * k h(x) are related to the coordinates z. More precisely, they correspond to coordinates ξ ξ ξ given by

ξ ξ ξ = Wz, W =     w * 1 . . . w * n     , (40)
i.e. they are linearizing coordinates, in which the dynamics is diagonal.

The nonlinear case

This section generalizes the connection between the eigenfunctions of Koopman operator and linear normal modes of vibration to NNMs as defined by Shaw and Pierre [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF]. Recall from Section 2 that NNM manifolds are characterized by the following properties:

1. the manifold is invariant; 2. it passes through a stable fixed point (the origin); 3. the plane tangent to it at the origin is a linear mode of the linearized dynamics.

Following the discussion of Section 4.1, consider the surface defined by the intersection of n -2 zero level sets of eigenfunctions

s k (x) = 0, ∀k = k 1 , k 2 , λ k 1 = λ k 2 . ( 41
)
Relying on the properties of the previous sections, we can show that the manifold defined by (41) verifies the properties 1-3:

1. As shown in Section 3.2.1 the level sets s k (x) = 0 are all invariant, therefore the intersection of any number of them is invariant as well. 2. From ( 35) and (37) it follows that

s k (0 0 0) = w * k h(0 0 0) = 0, (42)
thus any zero level set contains the origin, and so does their intersection. 3. The tangent space at the origin of (41) is defined by the linear equations

∂s k ∂x (0 0 0)x = 0, ∀k = k 1 , k 2 (43)
which, using (39), reduces to

w * k x = 0, ∀k = k 1 , k 2 , ( 44
)
and corresponds to the linear mode associated with λ k 1 and λ k 2 as shown in Section 4.1.

Therefore, the generalization of linear modes to nonlinear systems in the Koopman operator framework corresponds to the definition of NNMs proposed by Shaw and Pierre.

NNM approximation using Taylor series

Thanks to the connection established in the previous section, it is possible to use methods from the Koopman operator framework to compute NNMs. In this section a Taylor series approximation based on the spectral expansion ( 27) is proposed. The underlying idea is to use the values of the eigenfunctions as coordinates. This has the advantage that the dynamics in these coordinates is linear as illustrated in Section 4.2; furthermore, equations (41) have a particularly simple form in this coordinate system.

Let ξ k = s k (x) so that [START_REF] Neild | Nonlinear modal decomposition using normal form transformations[END_REF] becomes

id(x) = x = (k 1 ,...,kn)∈N n k 1 +...+kn>0 v k 1 ,...,kn ξ k 1 1 . . . ξ kn n , (45)
and ( 41)

ξ k = s k (x) = 0, ∀k = k 1 , k 2 , λ k 1 = λ k 2 . ( 46
)
Suppose k 1 = 1 and k 2 = 2. Equation (45) expresses the transformation between the coordinate system defined by the Koopman operator eigenfunctions ξ ξ ξ and the original coordinates x, while (46) defines the NNM in the ξ ξ ξ coordinates. Combining the two, we can express the points on the NNM as a function of the variables ξ 1 and ξ 2 :

x = (k 1 ,k 2 )∈N 2 k 1 +k 2 >0 v k 1 ,k 2 ξ k 1 1 ξ k 2 2 := Ψ Ψ Ψ(ξ 1 , ξ 2 ), (47)
where the notation v k 1 ,k 2 = v k 1 ,k 2 ,0,...,0 is used. Equation (47) assigns to any pair of values of the eigenfunctions s 1 (x) and s 2 (x) a point on the manifold, and each point on the manifold corresponds to a value ξ 1 = ξ 2 ∈ C. In other words, the manifold is parametrized using all the admissible values for the eigenfunctions. Moreover, it follows from ( 27) that the time evolution of a trajectory on the manifold is given by

x(t) = (k 1 ,k 2 )∈N 2 k 1 +k 2 >0 v k 1 ,k 2 ξ k 1 1 ξ k 2 2 exp((k 1 λ 1 + k 2 λ 2 )t) = Ψ Ψ Ψ(ξ 1 (t), ξ 2 t(t)). ( 48
)
Equation ( 47) can be truncated to a fixed order to approximate the invariant manifold in a neighborhood of the origin, and correspondingly (48) to approximate the NNM motion. The vectors v k 1 ,k 2 are the Koopman modes of the identity function and can be computed as discussed in Appendix A.

To conclude this section some intrinsic limitations of approximating NNMs with the above Taylor series are discussed. First of all, the eigenfunctions might be analytic only in a neighborhood of the origin, so that the amplitude to which NNMs can be computed using (47) might also be limited (see also [START_REF] Mauroy | Global stability analysis using the eigenfunctions of the Koopman operator, Automatic Control[END_REF] for a similar issue in approximating eigenfunctions with Taylor series). Even in cases in which the series converges in the whole basin of attraction, the computation of the Koopman modes is more demanding as the order increases, and becomes unfeasible for very high orders.

Due to this the method proposed in this section is not meant to be a ready-to-use algorithm for the computation of NNMs of arbitrary systems, but rather a simple method to compute NNMs for limited amplitudes in small-dimensional systems. Its main purpose is to show that the theory of Koopman operator can lead to computational methods that use the coordinates given by the eigenfunctions. More generic algorithms could be derived using the results of the following section, in which a system of PDEs resembling (3) is derived.

Uniqueness of NNMs

The section is concluded by a remark about conditions 1-3 of Section 4.3 defining the NNMs. Generally these conditions do not determine a unique manifold in the dissipative case. For the sake of illustration, consider a two-dimensional linear system

ẋ1 = σ 1 x 1 , ẋ2 = σ 2 x 2 , ( 56
)
with σ 2 < σ 1 < 0. The phase portrait, illustrated in Figure 3, shows that almost all trajectories approach the linear eigenspace x 2 = 0, leading to a continuum of tangent manifolds. Analogous phenomena are possible in higher dimensions. For completeness, Appendix B presents a family of invariant manifolds tangent to a linear normal mode, for a fourdimensional linear system. However, assuming nonresonant eigenvalues, the coefficients of a polynomial approximation are uniquely determined at any order, so that only one of the manifolds is analytic. This observation applies to Shaw and Pierre's approach as well (note that the existence of locally analytic manifolds is guaranteed by the existence of a local analytic change of coordinates, see Section 4.2). The issue of non-uniqueness is discussed in detail in [START_REF] Haller | Nonlinear normal modes and spectral submanifolds: Existence, uniqueness and use in model reduction[END_REF].

We note that also in the context of the Koopman operator there is the possibility to obtain multiple manifolds for only one eigenvalue. Specifically, our construction of the eigenfunctions relies on the analyticity of the transformation h of Section 4.2, which is uniquely defined by its Taylor series. However, if this smoothness condition is dropped it is possible to find different transformations verifying (35) and (36), obtaining in this way multiple eigenfunctions for the same eigenvalue.

A 2 DOF example

The Taylor series expansion derived in Section 4.4 and based on the Koopman operator framework is used to compute the NNMs of the system shown in Figure 2. The equations of motion are

m 1 ẍ1 + k a x 1 + c ẋ1 + k nl x 3 1 + k b (x 1 -x 2 ) + c( ẋ1 -ẋ2 ) = 0, m 2 ẍ2 + k a x 2 + c ẋ2 + k b (x 2 -x 1 ) + c( ẋ2 -ẋ1 ) = 0, ( 57
)
where m = m 2 = m = 1 Kg, c = 0.05 Ns/m, k a = 1 N/m and k nl = 0.5 N/m 3 . With

y 1 = ẋ1 , y 2 = ẋ2 , η 1 = (k a + k b )/m, η 2 = k b /m, β = c/m, α = k nl /m, (57) is recast in first-order form      ẋ1 ẋ2 ẏ1 ẏ2      =      0 0 1 0 0 0 0 1 -η 1 η 2 -2β β η 2 -η 1 β -2β           x 1 x 2 y 1 y 2      +      0 0 -αx 3 1 0      . ( 58
)
The Jacobian matrix

A =      0 0 1 0 0 0 0 1 -η 1 η 2 -2β β η 2 -η 1 β -2β      (59)
has two pairs of complex conjugate eigenvalues λ 1 = λ 2 , λ 3 = λ 4 , which verify a 3 : 1 resonance condition λ 1 = 3λ 3 for k b = 4. Exact linear resonance is eliminated by varying k b but folding is likely to occur in the neighborhood of those parameter values. Table 1 shows how the eigenvalues vary for different values of k b . To determine the coefficients describing the NNM in Eq. ( 47), the trajectory on the manifold given by ( 48) is injected into the equations of motion. Evaluating at time t = 0 yields

k b λ 1 f n,1 ζ 1 λ 2 f n,
(k 1 ,k 2 )∈N 2 k 1 +k 2 >0 v k 1 ,k 2 ξ k 1 1 ξ k 2 2 (k 1 λ 1 + k 2 λ 2 ) = A (k 1 ,k 2 )∈N 2 k 1 +k 2 >0 v k 1 ,k 2 ξ k 1 1 ξ k 2 2 +            0 0 -α         (k 1 ,k 2 )∈N 2 k 1 +k 2 >0 v k 1 ,k 2 ξ k 1 1 ξ k 2 2     1     3 0            , (60)
where (v) i stands for the ith of the vector v. Matching the terms corresponding to the first-order powers results in two systems

(λ 1 I -A)v 10 = 0, (λ 2 I -A)v 01 = 0, (61)
so that the first-order modes v 10 and v 01 are the eigenvectors of A corresponding to λ 1 and λ 2 . If there is no resonance between the eigenvalues of A, the successive systems admit a unique solution, and are obtained by matching like-power terms (see the appendix for the complete derivation). Doing so for the second-order terms yields linear equations of the form

((k 1 λ 1 + k 2 λ 2 )I -A)v k 1 k 2 = 0, k 1 + k 2 = 2, (62)
where the right hand side is zero due to the absence of second-order terms in the vector field. We can proceed iteratively, solving for any mode v k 1 ,k 2 a linear system of the form

(rI -A)v k 1 k 2 = b, r = k 1 λ 1 + k 2 λ 2 ( 63
)
where b is obtained as a sum of lower-order modes (see appendix A for the details). A change of coordinates using the Jordan basis is used to improve numerical stability.

The above procedure is used to compute the high-frequency (i.e. the out-of-phase) mode of the 2 DOF system for k b = 4.3. The Taylor series used to generate the figures of this section were truncated at the 50 th order. The maximum amplitude at which the manifolds are represented was chosen so that trajectories obtained with (48) agree with those computed through numerical integration. To have a quantitative measure of the difference between two trajectories the normalized mean square error (NMSE) was used:

NMSE = 100 σ 2 x K k=0 x(k) -x(k) 2 , ( 64
)
where x is the time series obtained with (48), x is the reference series obtained through numerical integration, and σ 2 x is the variance of the latter. The NMSEs obtained to validate the manifolds were always lower than 1%, a value usually taken as threshold for a good match between two time series.

The result is compared with the finite element (FE) method developed in [START_REF] Renson | An effective finite-element-based method for the computation of nonlinear normal modes of nonconservative systems[END_REF]. This numerical method is based on an explicit parametrization of the invariant surface using a pair of master coordinates and it is known to be accurate provided that a good mesh is employed. For the out-of-phase mode, the selected master coordinates correspond to the displacement and velocity of the first mass. As can be seen in Figure 4 the results from the two methods agree very well, but the FE method allows higher levels of amplitude to be reached, while the Taylor series validity is more limited.

The lower-frequency (in-phase) NNM was computed for the different values of k b reported in Table 1. The results are shown in Figure 5 and compared with the NNMs obtained using the FE algorithm. As k b is decreased from 4.7 (Figures 5a and5b) to 4.3 (Figures 5c and5d) and 4.1 (Figures 5e and5f), the manifold geometry changes and starts to fold, which prevents the use of master coordinates for the parametrization. This is also confirmed by Figure 6, where the manifold obtained for k b = 4.1 is shown with a viewpoint orthogonal to the x 2 y 2 plane. The FE algorithm reaches the folding but cannot go beyond, since it cannot converge in region where the master coordinates parametrization is not valid. In contrast the folding does not pose any problem with the parametrization obtained through the Koopman eigenfunctions.

In Figure 7 one of the trajectories considered for the validation is compared with the corresponding numerical integration, in the case k b = 4.1. As expected according to the values of the NMSE, the two trajectories cannot be distinguished by visual inspection.

Equation (48) also highlights a structure of the evolution that is reminiscent of that obtained using Fourier series with periodic trajectories. The main difference is that the number of terms in (48) is not the same at all orders, but increases with the order. Figure 8 shows the terms corresponding to the first three non-zero orders, for k b = 4.7 and k b = 4.1, and in Figure 9 their sum is compared with the trajectories obtained using numerical integration. When only terms up to order five are used the agreement between the two trajectories is not as good as in Figure 7, nonetheless, it is remarkable that the main aspects of the dynamics are retained by this approximation, at least in the case considered herein. It is also interesting to note that in the case k b = 4.7 the trajectory presents a lower harmonic content, and using the same order a better approximation is obtained.

Conclusion

This research investigates the intimate connection that exists between NNMs of damped systems and the Koopman operator associated with their dynamics. In particular, we have shown that there is a correspondence between the invariant manifolds defined by Shaw and Pierre and the zero level sets of particular eigenfunctions of the Koopman operator, and that this correspondence is a natural generalization of the linear case.

Through the Kooman operator eigenfunctions a new characterization of NNMs is possible, and a novel approach to their parametrization is proposed. The main advantage of this new approach is its global validity, obtained as a consequence of the linearizing properties provided by the eigenfunctions. This parametrization overcomes a known limitation in the use of master coordinates.

However, the Taylor series expansion used herein is local in nature and restricted to low-dimensional systems. In order to target arbitrarily high levels of amplitude, as well as high-dimensional systems, more efficient methods should be developed for solving the system of PDEs (51) or (55). As mentioned in Section 4.5 there are differences between these systems and the one obtained in the work of Shaw and Pierre. Nonetheless, all these systems share a common structure, which might allow the extension of numerical methods developed using master coordinates to this new framework. In addition, the Koopman operator approach is associated with recent data analysis techniques that are well-suited to high-dimensional systems [START_REF] Williams | A kernel approach to data-driven Koopman spectral analysis[END_REF]. This could potentially open the door to develop data-driven methods for computing NNMs of large systems.

The preliminary results presented here motivate a deeper investigation of NNMs for damped systems. For example it is well known that in conservative systems NNMs can bifurcate, and exceed the number of degrees of freedom of the system, but how these phenomena survive the presence of damping is still an open question.

|α| = n j=1 α j , (66)
∂ α f l = ∂ |α| f l ∂x α 1 1 . . . ∂x αn n , (67)
x α = n j=1 x α j j . (68)
Using this notation, we can write

s k 1 ,...,kn (x) = s k (x), k = (k 1 , . . . , k n ), (69)
with

s(x) =     s 1 (x) . . . s n (x)     . (70)
In addition, ( 26) can be rewritten as

x = h≥1 k∈N n |k|=h v k s k (x), (71) and
φ φ φ(t, x) t=0 = f (x) = h≥1 k∈N n |k|=h v k s k (x)kλ, (72)
where kλ = n j=1 k j λ j and λ j are the eigenvalues of the Jacobian matrix

A = ∂f ∂x x=0 0 0 . ( 73
)
Similarly, the Taylor series of f (x) can be written as

f l (x) = h≥1 k∈N n |k|=h 1 k! ∂ k f l (0 0 0)x k , (74)
and injecting (71) in (74) leads to

f l (x) = h 1 ≥1 k 1 ∈N n |k 1 |=h 1 1 k 1 ! ∂ k 1 f l (0 0 0) n j=1      h 2 >1 k 2 ∈N n |k 2 |=h 2 (v k 2 ) j s k 2 (x)      k 1 j . ( 75
)
By comparison of (72) and (75), for any k ∈ N n we obtain

(v k ) l kλ = |k| h=1 |j|=h 1 j! ∂ j f l (0 0 0) (k 1 ...k |j| )| i k i =k n i=1 i s js m= i-1 s js+1 (v km ) i . ( 76
)
These equations are linear in the terms with maximum value of |k|, and can be solved if all the v k with |k | < |k| are known. The matrix defining the linear system related to v k has the form kλI-A. The first-order modes are the eigenvectors of A and it is necessary to choose a scaling factor, which correspondingly fixes the value of the eigenfunctions (defined up to a multiplicative constant). If there are no resonances between the eigenvalues, the linear systems related to modes of higher order are uniquely determined, since they are defined by a nonsingular matrix. We can also use the known form of the matrices defining the linear systems to simplify the computation. In the case of a diagonalizable Jacobian matrix, the systems are diagonal when expressed in coordinates defined by a Jordan basis.

To obtain a manifold using (47) not all the modes are necessary, but just a subset of them: the v k , k = (k 1 , . . . , k n ), for which only two indices k j 1 and k j 2 are nonzero. Since to obtain the linear system corresponding to a given mode v k only the modes v k , k = (k 1 , . . . , k n ), with k l ≥ k l are needed, the computation can be restricted to the modes that appear in (47).

As concluding remarks we note that approximating a solution of (51) with a Taylor series yields the same linear systems. Also, the procedure presented here corresponds in essence to a normal form approximation. In that context, starting from a dynamical system with diagonal linear part, a transformation can be obtained so that the resulting equations of motion are the simplest possible. The only algebraic condition needed to obtain a transformation that linearizes the dynamic is the nonresonant eigenvalues. The corresponding linearizing coordinates are the eigenfunctions of the Koopman operator.

B. A family of invariant manifolds for a linear system

Consider a four-dimensional linear system ẋ = Ax,

with

A =      -σ 1 ω 1 0 0 -ω 1 -σ 1 0 0 0 0 -σ 2 ω 2 0 0 -ω 2 -σ 2      . ( 78
)
Note that any four-dimensional linear system with two couples of complex conjugate eigenvalues can be transformed into this form. It will be assumed that σ 1 > 0 and σ 2 > 0 (stability).

For any value of c 1 and c 2 consider the manifold x(u, v) = f (u, v; c 1 , c 2 ),

f 1 (u, v) = u f 2 (u, v) = v f 3 (u, v) = u 2 + v 2 σ 2 2σ 1 cos log (u 2 + v 2 ) ω 2 2σ 1 c 1 -sin log (u 2 + v 2 ) ω 2 2σ 1 c 2 f 4 (u, v) = u 2 + v 2 σ 2 2σ 1 sin log (u 2 + v 2 ) ω 2 2σ 1 c 1 + cos log (u 2 + v 2 ) ω 2 2σ 1 c 2 . ( 79
)
The derivatives of f are

f u,1 = 1 f u,2 = 0 f u,3 = u (u 2 + v 2 ) σ 2 2σ 1
-1 cos log(u 2 +v 2 )ω2

σ 1 f v,1 = 0 f v,2 = 1 f v,3 = v (u 2 + v 2 ) σ 2 2σ 1 -1 cos log(u 2 +v 2 )ω2 2σ 1 (c 1 σ 2 -c 2 ω 2 ) -sin log(u 2 +v 2 )ω2 2σ 1 (c 2 σ 2 + c 1 ω 2 ) σ 1 f v,4 = v (u 2 + v 2 ) σ 2 2σ 1 -1 cos log(u 2 +v 2 )ω2 2σ 1 (c 2 σ 2 + c 1 ω 2 ) + sin log(u 2 +v 2 )ω2 2σ 1 (c 1 σ 2 -c 2 ω 2 ) σ 1 . (80)
These are continuous for σ 2 > σ 1 , and more regular if the distance between the two real parts of the eigenvalues is increased further. Assuming sufficient regularity, f defines a manifold, with tangent plane at a point (u, v) spanned by the vectors

f u = ∂f ∂u , f v = ∂f ∂v . ( 81
)
Since it holds Af (u, v) = (-σ

1 u + ω 1 v)f u + (-ω 1 u -σ 1 v)f v , (82)
all the manifolds are invariant ((82) is the analogous of the equation derived by Shaw and Pierre equation for the system (77)). Moreover, if σ 2 > σ 1

f u (0, 0) =      1 0 0 0      , f v (0, 0) =      0 1 0 0      , ( 83
)
so that at the origin the tangent plane is x 3 = 0, x 4 = 0, which corresponds to the linear normal mode, but only for c 1 = c 2 = 0 the whole manifold is the linear normal mode. Note that apart from the points verifying x 1,0 = x 2,0 = 0 (i.e. points on the fast linear mode), for any other point (x 1,0 , x 2,0 , x 3,0 , x 4,0 ), there exists a manifold which passes through it, this is obtained using the parameters (84)

c 1 = x 3,0 1 x 2 1,0 + x 2
As a final remark we note that even if each couple of coefficients c 1 , c 2 defines a different manifold, the derivatives of all the functions defined by (79) have the same values at the origin (up to the order at which they are defined): they all share the same polynomial approximation. Interestingly the same situation is encountered in the context of center manifolds.

2 Figure 3 :

 23 Figure 3: Phase portrait of a two-dimensional linear system

Figure 4 :

 4 Figure 4: Out-of-phase NNM; FE method, Taylor series; (a) x 1 , y 1 , x 2 ; (b) x 1 , y 1 , y 2

Figure 5 : 2 Figure 6 :

 526 Figure 5: In-phase NNM for different values of k b ; (a)-(b) k b = 4.7, (c)-(d) k b = 4.5, (e)-(f) k b = 4.1; FE method, Taylor series 19

Figure 7 :

 7 Figure 7: Comparison of time series; Numerical integration, Taylor approximation; initial condition: x 1 (0) = -0.23, x 2 (0) = -0.76, y 1 (0) = -0.4, y 2 (0) = -0.86; (a) x 1 , (b) x 2 , (c) y 1 , (d) y 2

Figure 8 :

 8 Figure 8: Decomposition of trajectories in different terms; 1 st order term, 3 rd order term, 5 th order term; (a) k b = 4.7 (x 1 (0) = -0.57, x 2 (0) = -0.55, y 1 (0) = -0.44, y 2 (0) = -0.66); (b) k b = 4.1 (x 1 (0) = -0.23, x 2 (0) = -0.76, y 1 (0) = -0.4, y 2 (0) = -0.86)

Figure 9 :

 9 Figure 9: Comparison of 5 th order approximation and numerical integration; Numerical integration, 5 th order approximation; (a) k b = 4.7 (x 1 (0) = -0.57, x 2 (0) = -0.55, y 1 (0) = -0.44, y 2 (0) = -0.66), (b) k b = 4.1 (x 1 (0) = -0.23, x 2 (0) = -0.76, y 1 (0) = -0.4, y 2 (0) = -0.86)

Table 1 :

 1 Eigenvalues, natural frequencies and damping ratios for different values of k b .

	2	ζ 2

  (c 1 σ 2 -c 2 ω 2 ) -sin log(u 2 +v 2 )ω2 (c 2 σ 2 + c 1 ω 2 ) + sin log(u 2 +v 2 )ω2 2σ 1 (c 1 σ 2 -c 2 ω 2 )

					2σ 1	2σ 1	(c 2 σ 2 + c 1 ω 2 )
					σ 1
	f u,4 =	u (u 2 + v 2 )	σ 2 2σ 1	-1 cos	log(u 2 +v 2 )ω2 2σ 1

Acknowledgment

This paper presents research results of the Belgian Network DYSCO (Dynamical Systems, Control, and Optimization), funded by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office.

This work was performed while A. Mauroy held a return grant from the Belgian Science Policy (BELSPO). The author L. Renson is a Marie-Curie COFUND Postdoctoral Fellow of the University of Liège, co-funded by the European Union.

Analogy with Shaw and Pierre's manifold governing PDE

In Section 2 a system of PDEs (3) is derived starting from Shaw and Pierre definition of NNMs. A similar system can be derived for Ψ Ψ Ψ(ξ 1 , ξ 2 ) defined in (47). To this end, note that the dynamics of ξ k is by definition ξk

It follows that the time derivative of Ψ Ψ Ψ(ξ 1 (t), ξ 2 (t)) is given by

Substituting (50) into the equations of motion [START_REF] Shaw | Normal modes for non-linear vibratory systems[END_REF], and evaluating at t = 0 yields

Equation ( 51) is a set of n PDEs very similar in principle to the one derived in Section 2. The approaches followed here and in Section 2 can be both interpreted as a way to solve

When using state-space variables as in (3), Eq. ( 52) is inverted, but the inversion might be valid only locally. In contrast, a different set of coordinates is used in (51) to parametrize the solution, which avoids issues arising when NNMs have a complex geometry (e.g. folding).

The two systems of PDEs (3) and (52) have a similar structure, but (3) corresponds to n -2 equations while (52) consists of n equalities. Moreover, the second system admits more solutions, due to the fact that the eigenfunctions are defined modulo a scaling factor.

) is a solution for any complex c = 0.

To obtain Ψ Ψ Ψ(ξ 1 , ξ 2 ) real-valued ξ 1 = ξ 2 must be assumed in (51). Equivalently, a system of equations can be derived using as variables u = Re(ξ 1 ) and v = Im(ξ 1 ), whose dynamics is

with

with ψ ψ ψ(u, v) = Ψ Ψ Ψ(u + iv, uiv). Clearly, the above procedure can be followed for any pair of complex conjugate eigenfunctions, yielding a system of PDEs (51) or (55) for each NNM. These systems differ only for the value of λ (or σ and ω).

A. Computation of the Koopman Modes of the Identity

In this appendix we derive a method for computing the Koopman modes v k 1 ,...,kn of the identity (see [START_REF] Neild | Nonlinear modal decomposition using normal form transformations[END_REF]). In particular, we show that these modes are the solutions of linear systems.

We fix some notation. For α = (α 1 , . . .


Enke Baldini

Chiara Tuccilli

Yannick Arlot-Bonnemains

Frank Chesnel

Salvatore Sorrenti

Corrado De

Antonio Catania

Eleonora D ' Armiento

Alessandro Antonelli

Poupak Fallahi

Sara Watutantrige-Fernando

Francesco Tartaglia

Susi Barollo

Caterina Mian

Marco Bononi

Stefano Arceri

Domenico Mascagni

Massimo Vergine

Daniele Pironi

Massimo Monti

Angelo Filippini

Salvatore Ulisse
email: salvatore.ulisse@uniroma1.it

Prof Salvatore Ulisse

Deregulated expression of VHL mRNA variants in papillary thyroid cancer

Keywords: Papillary thyroid cancer, VHL, Gene expression, Prognosis, Biomarkers

Recent findings demonstrated that a subset of papillary thyroid cancers (PTCs) is characterized by reduced expression of the von Hippel-Lindau (VHL) tumor suppressor gene, and that lowest levels associated with more aggressive PTCs. In the present study, the levels of the two VHL mRNA splicing variants, VHL-213 (V1) and VHL-172 (V2), were measured in a series of 96 PTC and corresponding normal matched tissues by means of quantitative RT-PCR. Variations in the mRNA levels were correlated with patients' clinicopathological parameters and disease-free interval (DFI).

The analysis of VHL mRNA in tumor tissues, compared to normal matched tissues, revealed that its expression was either up-or down-regulated in the majority of PTC. In particular, V1 and V2 mRNA levels were altered, respectively, in 78 (81.3%) and 65 (67.7%) out of the 96 PTCs analyzed. A significant positive correlation between the two mRNA variants was observed (p<0.001). Univariate analysis documented the lack of association between each variant and clinicopathological parameters such as age, tumor size, histology, TNM stage, lymph node metastases, and BRAF mutational status. However, a strong correlation was found between altered V1 or V2 mRNA levels and DFI. Multivariate regression analysis indicated higher V1 mRNA values, along with lymph node metastases at diagnosis, as independent prognostic factors predicting DFI. In conclusion, the data reported demonstrate that VHL gene expression is deregulated in the majority of PTC tissues. Of particular interest is the apparent protective role exerted by VHL transcripts against PTC recurrences.

Introduction

The von Hippel-Lindau (VHL) gene, located on chromosome 3p25-p26, consists of 3 exons [START_REF] Latif | Identification of the von Hippel-Lindau disease tumor suppressor gene[END_REF]. Two different VHL mRNAs are generated because of differential splicing: the variant 1 (V1), including all three exons, and the variant 2 (V2), lacking exon 2 (Iliopolus et. 1998;[START_REF] Gnarra | Mutations of the VHL tumour suppressor gene in renal carcinoma[END_REF][START_REF] Richards | Expression of the von Hippel-Lindau disease tumour suppressor gene during human embryogenesis[END_REF]. From the V1 mRNA, two VHL protein (pVHL) isoforms are synthesized due to alternate initiation sites of translation; the pVHL 213 , of about 30 kDa and containing 213 amino acids, and the pVHL 160 , of about 19 kDa and containing 160 amino acids [START_REF] Schoenfeld | A second major native von Hippel-Lindau gene product, initiated from an internal translation start site, functions as a tumor suppressor[END_REF]. The V2 mRNA variant encodes a protein of 172 amino acids lacking an amino-terminal pentameric acid repeat (β-domain), with a predicted molecular weight of about 19 kDa (pVHL 172 ) [START_REF] Chesnel | The von Hippel-Lindau tumour suppressor gene: uncovering the expression of the pVHL172 isoform[END_REF].

The VHL gene is widely expressed in human tissues and acts as a tumor suppressor gene [START_REF] Los | Expression pattern of the von Hippel-Lindau protein in human tissue[END_REF]. Loss-of-function mutations cause the so-called VHL disease, an autosomal dominant disorder characterized by retinal angioma, cerebellar and spinal hemangioblastomas, clear-cell renal cell carcinoma (ccRCC), pheochromocytoma and pancreatic neuroendocrine tumor [START_REF] Maher | von hippel-Lindau disease: a clinical and scientific review[END_REF]. Besides, a number of sporadic cancers, including ccRCC, are strongly associated with VHL reduced expression and/or loss-of-function mutations [START_REF] Gnarra | Mutations of the VHL tumour suppressor gene in renal carcinoma[END_REF][START_REF] Cassol | Endocrine manifestations of von Hippel-Lindau disease[END_REF]. The best characterized pVHL function is ubiquitination followed by proteasome degradation of target proteins [START_REF] Robinson | The multifaceted von Hippel-Lindau tumor suppressor protein[END_REF]. In particular, the pVHL forms a multiprotein complex with elongins B and C, cullin 2 and Rbx-1, which functions as a ubiquitinligating enzyme (E3 ligase). Within the complex pVHL is responsible for recognition of substrates, among which the members of the hypoxia-inducible factor α (HIFα) family [START_REF] Robinson | The multifaceted von Hippel-Lindau tumor suppressor protein[END_REF]. In complex with HIFβ, HIFα act as transcription factors enhancing the expression of a variety of genes involved in the adaptive response to hypoxic condition and tumor progression as well, such as genes that promote neoangiogenesis (e.g. vascular endothelial growth factor, VEGF), energy metabolism (e.g. glucose transporter 1, GLUT1), erythropoiesis (e.g. erythropoietin, EPO), and cell survival (e.g. transforming growth factor-α, TGF-α) [START_REF] Robinson | The multifaceted von Hippel-Lindau tumor suppressor protein[END_REF][START_REF] Balamurugan | HIF-1 at the crossroads of hypoxia, inflammation, and cancer[END_REF]. In addition, VHL has been shown to affect several hypoxia-independent cellular functions, including extracellular matrix formation, spindle microtubule stability, cilia formation, epithelial-to-mesenchymal (EMT) transition, cell proliferation, apoptosis and DNA damage response [START_REF] Robinson | The multifaceted von Hippel-Lindau tumor suppressor protein[END_REF].

The pVHL region interacting with the elongins and Cul-2, named BC box, is coded by the third exon, thus all pVHL isoforms are able to form the E3 ligase complex (UniProtKB -P40337).

However, the binding with HIFα occurs via the β-domain, which is missing in the pVHL 172 [START_REF] Bonicalzi | Role of exon 2-encoded beta -domain of the von Hippel-Lindau tumor suppressor protein[END_REF]. Such domain is also essential for VHL interaction with the chaperonincontaining t-complex polypeptide 1 (CCT), a cytosolic molecular chaperone that assists in the folding of actin, tubulin and other cytosolic proteins. Even if the physiological role of pVHL 172 remains to be elucidated, it might be supposed that this isoform behaves, at least partially, as a VHL 213 antagonist, and possibly is endowed with new functions. Some reports suggested that pVHL could play a role also in the progression of epithelial thyroid cancer (TC) (The Cancer Genome Atlas Research Network, 2014;[START_REF] Stanojevic | Low VHL mRNA expression is associated with more aggressive tumor features of papillary thyroid cancer[END_REF][START_REF] Hinze | Expression of the von Hippel-Lindau tumor suppressor gene in nonneoplastic and neoplastic lesions of the thyroid[END_REF][START_REF] Hunt | Loss of heterozygosity of the VHL gene identifies malignancy and predicts death in follicular thyroid tumors[END_REF]. The latter represents the most common endocrine malignancy accounting for roughly 1% of all human cancers, and its incidence has been increasing over the last decades mainly due to the improved ability to diagnose malignant transformation in small nonpalpable thyroid nodules [START_REF] Jemal | Cancer Statistics, 2009[END_REF][START_REF] Davies | Increasing incidence of thyroid cancer in the United States, 1973-2002[END_REF][START_REF] Kinder | Well differentiated thyroid cancer[END_REF][START_REF] Patel | Poorly differentiated and anaplastic thyroid cancer[END_REF][START_REF] Pasieka | Anaplastic thyroid cancer[END_REF]. More than 90% of TC are differentiated thyroid carcinomas (DTC), 70-80% of which are represented by the papillary (PTC), and the remaining by the follicular (FTC) histotype [START_REF] Nikiforov | Diagnostic pathology and molecular genetics of the thyroid[END_REF]. Although derived from the same cell type, the DTC show specific histological features, biological behavior and degree of differentiation because of different genetic alterations (The Cancer Genome Atlas Research Network, 2014;[START_REF] Nikiforov | Diagnostic pathology and molecular genetics of the thyroid[END_REF]). Among the somatic activating mutations, those of genes involved in the mitogen activated protein kinase (MAPK) signaling pathway, e.g. Ras, BRAF and RET/PTC rearrangements, are held responsible for the majority of PTC (The Cancer Genome Atlas Research Network, 2014;[START_REF] Nikiforov | Diagnostic pathology and molecular genetics of the thyroid[END_REF]. Treatment of patients includes total thyroidectomy followed, if necessary, by 131 I therapy.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
The prognosis is generally favorable, with a 10-year-survival rate of approximately 90%.

Nevertheless, nearby 20% of patients have disease recurrences and tumor-related deaths [START_REF] Eustatia-Rutten | Survival and death causes in differentiated thyroid carcinoma[END_REF]. The stratification and prognosis of patients depends on clinicopathological variables such as age, tumor size, histology, lymph nodal or distant metastases [START_REF] Eustatia-Rutten | Survival and death causes in differentiated thyroid carcinoma[END_REF][START_REF] Gospodarowicz | Prognostic factors in cancer[END_REF][START_REF] Passler | Prognostic factors of papillary and follicular thyroid cancers: differences in an iodine-replete endemic goiter region[END_REF][START_REF] Castagna | Delayed risk stratification, to include the response to initial treatment (surgery and radioiodine ablation), has better outcome predictivity in differentiated thyroid cancer patients[END_REF]). These parameters, however, are capable of providing only a rough prediction of the disease outcome placing patients with very different disease-specific progression and survival times within the same risk group.

Similarly, they fail to predict the risk of cancer relapse. Therefore, the identification of new prognostic molecular biomarkers able to testify tumor aggressiveness is required [START_REF] Handkiewicz-Junak | Molecular prognostic markers in papillary thyroid cancer: current status and future directions[END_REF][START_REF] Baldini | The urokinase plasminogen activating system in thyroid cancer: clinical implications[END_REF][START_REF] Ulisse | High expression of the urokinase plasminogen activator and its cognate receptor associates with advanced stages and reduced disease-free interval in papillary thyroid carcinoma[END_REF].

Although, to date, no loss-of-function mutations of the VHL gene in PTC tissues have been described, recent studies demonstrated that a subset of PTCs was characterized by low levels of VHL mRNA, which associated with more aggressive PTCs (The Cancer Genome Atlas Research Network, 2014;[START_REF] Stanojevic | Low VHL mRNA expression is associated with more aggressive tumor features of papillary thyroid cancer[END_REF].

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

In the present work, we measured the expression of the VHL gene, at mRNA level, in a case-study of 96 PTC tissues compared with their normal matched counterparts, and we evaluated the correlation of the VHL expression changes with clinicopathological parameters and disease-free interval.

Patients and Methods

Tissue samples, histology and patient's staging

Normal and matched tumor thyroid tissues were obtained from surgical specimens of 96 patients substitution at nucleotide 1799 (V600E), using the procedure previously described [START_REF] Barollo | BRAF in primary and recurrent papillary thyroid cancer: the relationship with (131)I and 2[END_REF].

Extraction of mRNA and quantitative RT-PCR

Frozen normal and tumor thyroid tissues were homogenized with the ultra-turrax, and total RNA extracted applying the acid guanidinium thiocyanate-phenol-chloroform method [START_REF] Chomczynsky | Single step method of RNA isolation by guanidinium thiocyanate-phenol-chloroform extraction[END_REF]. The first cDNA strand was synthesized from 5 µg of RNA with M-MLV reverse transcriptase and anchored oligo(dT)23 primers (Sigma Chemicals Co.). Parallel controls for DNA contamination were carried out omitting the reverse transcriptase. The templates obtained were used for quantitative PCR amplifications of the VHL mRNA variants 1 and 2, and three different housekeeping genes (GAPDH, RPL13A and SDHA) employing the LightCycler instrument (Roche Diagnostics, Mannheim, Germany), the SYBR Premix Ex Taq II (TliRNase H Plus) (Takara, Otsu, Shiga, Japan) and specific primers listed in Table 1.

Table 1. Sequences, genomic positions, and amplicon sizes of the primers used in qRT-PCR for the target and reference genes. VHL, von Hippel-Lindau; GAPDH, glyceraldehyde-3phosphate dehydrogenase; RPL13A, ribosomal protein L13a; SDHA, succinate dehydrogenase complex, subunit A; V1, VHL variant 1; V2, VHL variant 2.

Gene

Statistical analysis

First, the Shapiro-Wilk test was used to check whether the mRNA data were normally distributed, and they were not. Thus, the non-parametric Mann-Whitney U-test was used to calculate the statistical significance of differences in the expression levels of VHL mRNA V1 and V2 in female

vs

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

For the latter, values were classified based on the following criteria: fold change > 1.2 as "increased"; fold change < 0.8 as "decreased"; 0.8 ≤ fold change ≤ 1.2 as unvaried. All the statistical analyses were carried out with the SPSS software (IBM, Armonk, NY, USA), and the results were considered significantly different if p values were lower than 0.05.

Results

VHL gene expression in papillary thyroid cancer (PTC) tissues

The analyses of V1 and V2 mRNA levels in PTC tissues, compared to their normal matched tissues, revealed that V1 mRNA was deregulated in 78 out of 96 PTC (81.3%), with an increase in 36 and a decrease in 42 cases (Figure 1A). The V2 mRNA was altered in 65 out of 96 PTC (67.7%), being up-regulated in 24 and down-regulated in 41 cases (Figure 1B). As reported in Figure 1C, mRNA levels of the two VHL variants positively correlated to each other (p<0.001). Absolute quantification of the V1 and V2 mRNA copy number for each tumor and normal matched tissue showed that V1 mRNA copy number was higher compared to the V2 mRNA, with a median V1/V2 ratio of 3.95 in normal tissues (Figure 2). This ratio did not change significantly in PTC tissues, the median being 3.52 (p=0.655).

BRAF V600E mutation and VHL expression in PTC tissues

To assess the effect of BRAF V600E mutation on VHL expression, we analyzed the V1 and V2 mRNA levels in 76 PTC tissues for which the analysis of BRAF gene status was available. Of these, 37 PTC (48.7%) harbored the BRAF V600E mutation while 39 (51.3%) had the wild type BRAF. The results, reported in table 2, showed that the BRAF V600E mutation did not affect VHL gene transcription in PTC tissues.

Prognostic relevance of VHL expression in PTC patients.

Among the clinicopathological parameters analyzed, a significant association emerged between tumor relapse and low levels of V2, and a similar trend (p=0.057), was recorded for V1 (Table 2).

Moreover, increased V1 mRNA levels were found in male PTC tissues compared to the female ones (p=0.03). Cox regression analysis indicated that, among all the parameters considered, the presence of lymph node metastases at diagnosis represented a risk factor for tumor recurrences (Hazard Ratio, HR 9.37, p=0.003), while patients with increased V1 levels had a reduced risk of tumor recurrences (HR 0.18,p=0.025). Kaplan-Mayer analysis demonstrated a significant correlation of both V1 and V2 mRNA levels with patients' disease-free interval (DFI). In figure 3, panels A and

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

Discussion

To date, relative few studies have attempted to investigate the effects of a deregulated VHL expression in DTC evolution (The Cancer Genome Atlas Research Network, 2014;[START_REF] Stanojevic | Low VHL mRNA expression is associated with more aggressive tumor features of papillary thyroid cancer[END_REF][START_REF] Hinze | Expression of the von Hippel-Lindau tumor suppressor gene in nonneoplastic and neoplastic lesions of the thyroid[END_REF][START_REF] Hunt | Loss of heterozygosity of the VHL gene identifies malignancy and predicts death in follicular thyroid tumors[END_REF]. Earlier immunohistochemical analyses demonstrated the presence of pVHL in normal thyroid follicular epithelium as well as in non-neoplastic lesions, and in DTC [START_REF] Hinze | Expression of the von Hippel-Lindau tumor suppressor gene in nonneoplastic and neoplastic lesions of the thyroid[END_REF]. Intriguingly, an inverse correlation between VHL expression and tumor differentiation was found, the pVHL being diminished in poorly differentiated TC, and very weak or undetectable in the most deadly type of TC, the anaplastic carcinoma (Hinze et al. 2003).

However, no correlation between pVHL and VEGF proteins could be observed, suggesting a HIFαindependent role of pVHL in TC progression (Hinze et al. 2003). Although no loss-of-function mutations of VHL have been described so far in PTC, one study reported loss of heterozygosis (LOH) at the VHL gene in FTC [START_REF] Hunt | Loss of heterozygosity of the VHL gene identifies malignancy and predicts death in follicular thyroid tumors[END_REF]. These authors identified LOH in malignant but not in benign follicular tumors, and described a significant association between LOH and disease recurrence [START_REF] Hunt | Loss of heterozygosity of the VHL gene identifies malignancy and predicts death in follicular thyroid tumors[END_REF]. In PTC, no LOH specific for the VHL gene has been identified so far; however, different studies documented a variable LOH for the chromosomal region 3p where the VHL gene is located [START_REF] Grebe | Frequent loss of heterozygosity on chromosomes 3p and 17p without VHL or p53 mutations suggests involvement of unidentified tumor suppressor genes in follicular thyroid carcinoma[END_REF][START_REF] Rodrigues-Serpa | Loss of heterozygosity in follicular and papillary thyroid carcinomas[END_REF]. In particular, Grebe and colleagues reported the presence of LOH on the chromosome 3p in 28.5% of PTC analyzed, while Rodrigues-Serpa made the same observations in 40% of PTC [START_REF] Grebe | Frequent loss of heterozygosity on chromosomes 3p and 17p without VHL or p53 mutations suggests involvement of unidentified tumor suppressor genes in follicular thyroid carcinoma[END_REF][START_REF] Rodrigues-Serpa | Loss of heterozygosity in follicular and papillary thyroid carcinomas[END_REF]. Since in our PTC series we found that both V1 and V2 mRNAs were downregulated in about 43% of cases, it may be speculated that haploinsufficiency is, at least in part, the cause of this reduction. Nevertheless, further studies will be required to prove this point.

More recently, Stanojevic and colleagues (2015) evaluated the status and expression of the VHL gene on a case study of 264 patients with PTC. No somatic mutations or evidence of VHL downregulation via promoter hypermethylation were found. However, low VHL mRNA levels strongly associated with older patient's age, advanced clinical stage, classical PTC variant, and multifocality [START_REF] Stanojevic | Low VHL mRNA expression is associated with more aggressive tumor features of papillary thyroid cancer[END_REF]. The authors noticed also a marginal influence of low VHL expression on disease-free interval (DFI) (p = 0.0502).

Over the years, most of the studies regarding VHL expression have been accomplished without distinguishing between isoforms. However, as mentioned above, the absence of part of the βdomain (aa 114-154) in pVHL 172 modifies the number of beta sheets in the structure, which is likely to cause altered protein activity and interactions [START_REF] Schoenfeld | A second major native von Hippel-Lindau gene product, initiated from an internal translation start site, functions as a tumor suppressor[END_REF]. In our case study we sought to examine separately VHL isoforms in order to determine any differences in their expression profile that could be of relevance in thyroid tumor progression. Unfortunately, this approach implies the impossibility to discriminate isoforms at the protein level by means of showed, however, the lack of any influence of the BRAF V600E mutation on VHL expression.

A significant association of low V2 mRNA levels with tumor recurrence was found, and a similar trend (p=0.057) appeared for the V1 variant. In agreement with this, multivariate analysis demonstrated that patients bearing increased levels of V1 mRNA had a reduced risk of tumor recurrences. In addition, the output of Kaplan-Meier analysis evidenced a shorter DFI for patients with reduced level of the V2 mRNA in tumor tissues, and for patients with either unvaried or diminished V1 mRNAs. The protective effect of the VHL gene in PTC patients is in line with its anti-oncogenic role extensively documented for kidney cancer.

The discrepancy between our observations and those reported by Stanojevic and colleagues (2015) may be explained by the fact that they analyzed VHL transcripts without distinguishing the two splicing variants, and they adopted a different method of relative quantification based on the comparison of PTC samples with each other, rather than with the normal matched samples.

Anyway, while taking into account the conflicting statistical results of the two studies, Stanojevic and colleagues pointed to an inverse correlation between VHL expression and PTC aggressive behaviour, essentially reflecting the overall indications emerged in our findings.

In conclusion, although the data here reported need to be confirmed by means of larger case-studies, they demonstrate that the expression of both VHL transcriptional variants is deregulated in the majority of PTC tissues. Their role(s) in PTC progression remains, however, to be clarified.

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

ACCEPTED MANUSCRIPT

(

  19 males and 77 females, age range 11-83 yr, median 44 yr) who underwent total thyroidectomy for PTC at the Department of Surgical Sciences, "Sapienza" University of Rome (39 patients) or at the Department of Medicine, University of Padua (57 patients). All patients gave their informed consent, and the study was approved by the local ethical committee (Protocol No. 2615). Tissue samples were collected, frozen in liquid nitrogen and stored at -80 °C. Of the 96 PTC patients, 72 exhibited classical, 17 follicular, 3 tall cell and 4 oncocytic variants. The histological diagnoses were made independently by two different histopathologists according to the World Health Organization classification (Hedinger et al. 1989). At the time of surgery lymph node metastases were found in 39 patients. Following TNM staging, 62 patients were at stage I, 1 at stage II, 27 at stage III and 6 at stage IV. Approximately 40 to 50 days later, all the patients underwent radioiodine treatment followed by thyroid hormone replacement therapy. To ascertain their disease-free condition, 4 to 5 months later all the patients underwent neck ultrasound and serum Tg measurement. Recurrences were diagnosed by measurement of serum Tg levels either in basal conditions or following recombinant human TSH stimulation; FNA cytology and/or Tg determination in the FNA wash-out from lymph nodes; 131 I whole body scan; histological analysis following surgical resection of the lesion. The follow-up included 80 patients (mean 57.9±36.0 months, range 5-141 months), 54 of whom were at TNM stage I. During the follow-up 17 recurrences were recorded. Determination of BRAF V600E mutation Genomic DNA was extracted from the frozen tumor tissues using the DNeasy Blood and Tissues kit (QIAGEN, Milan, Italy) according to the manufacturer's protocol. The BRAF status of exon 15 was assessed by both direct sequencing and mutant allele-specific PCR amplification for the T to A M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT 6

  checked by automated DNA sequencing (Bio-Fab Research, Rome, Italy), evaluation of melting temperatures, and electrophoresis on 2% agarose gel containing ethidium bromide. For the relative quantification of gene expression, standard curves were made with five-fold dilutions of mixed cDNAs from human thyroid tissues for all reference and target genes. Calculation of data was performed by the Relative Expression Software Tool (REST 2009) using a normalization factor computed as geometric media of the 3 reference genes, as previously described[START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Ulisse | In papillary thyroid carcinoma BRAFV600E is associated with increased expression of the urokinase plasminogen activator and its cognate receptor, but not with disease-free interval[END_REF]. The fold change of VHL gene expression for each tumor sample was referred to its normal counterpart.For the absolute quantification of VHL V1 and V2 mRNAs, two standard curves were created with five-fold dilutions of known amounts of two different BamHI-linearized plasmids, each containing the sequence of a VHL variant. Specifically, pCDNA3.1-FlagHA-VHL213 or -VHL172 plasmids were generated by inserting a 639bp or a 516bp BamHI-XhoI fragment (coding for human VHL open reading frames digested out of pET21-VHL213 or pET21-VHL172 plasmids, respectively)[START_REF] Chesnel | The von Hippel-Lindau tumour suppressor gene: uncovering the expression of the pVHL172 isoform[END_REF] in frame to the 3'-end of the Flag-HA sequence of pcDNA3-FlagHA plasmid (a kind gift from Dr S. Rouquier, Toulouse, France). Constructs were confirmed by Sanger sequencing, amplified in NEB5alpha bacteria (New England Biolabs, Evry, France) and purified using the Nucleobond Endofree plasmid purification kit (Macherey Nagel, Hoerdt, France). Copy numbers of the VHL mRNA variants were estimated by interpolation of the crossing points obtained for the tissue samples on these standard curves, and quantitative ratios VHL V1/V2 were calculated for normal and tumor tissue of each patient.

Figure 1 :

 1 Figure 1: Relative expression levels of the VHL mRNA splicing variants in 96 PTC tissues. A) and B) Variations of the V1 and V2 mRNAs, respectively. The fold changes were calculated considering the mRNA values of normal matched thyroid tissues equal to one. The statistical evaluation of the data was performed with the non-parametric Mann-Whitney test. The small bars in the graph indicate the median values. C) Correlation analysis of V1 and V2 mRNAs in PTC tissues. The data were evaluated by applying the Rho Spearman test.

Figure 2 :

 2 Figure 2: VHL V1/V2 mRNAs ratio in PTC and normal matched tissues.

  DFI is represented for patients grouped in three categories based on the mRNA fold changes: increased, unvaried and decreased. Looking at the graphic of V1 mRNA, the trends for patients with unvaried and decreased values were very similar, so they were grouped together and compared against PTC patients with increased values. As shown in figure3 (panel C), the returned p-value of log-rank test was 0.007. Regarding the V2 mRNA, PTC patients with unvaried or increased values had analogous DFI profiles (figure3, panel B). When the latter were grouped together and compared to patients with decreased V2 mRNA levels, a significant difference (p=0.025

Figure 3 :

 3 Figure 3: Levels of VHL mRNA variants and disease-free interval (DFI) in PTC patients. Kaplan-Mayer analysis combined with Mantel-Cox log-rank statistical test, performed on 80 PTC patients followed-up from 5 to 141 months. Expression values of V1 and V2 VHL variants were classified as increased (fold change > 1.2), decreased (fold change < 0.8) or unvaried (0.8 ≤ fold change ≤ 1.2).

  male patients; in classical PTC variant vs other variants; in BRAF V600E mutated vs wild type

	PTC; in metastatic (N1) vs non-metastatic (N0) PTC; in T1-2 vs T3-4 tumor sizes; in TNM I-II vs
	III-IV stages; in presence or absence of recurrence. The correlation between V1 and V2 mRNAs,
	and between each of them and patient's age was evaluated by the Spearman's Rho test.
	The VHL V1/V2 mRNA ratio comparison between normal and PTC tissues was performed with the
	Wilcoxon signed-rank test. To assess the independent association of patient's age, gender, tumor

size, histological variants, BRAF status, lymph node metastases, stage and VHL mRNA variants with recurrences, the Cox regression with backward analysis was used. The impact of V1 or V2 expression on DFI was assessed by the Kaplan-Meier analysis combined with Mantel-Cox log-rank.

Table 2 . Univariate statistical analysis of VHL mRNA variants, and PTC patient's characteristics and high-risk clinicopathological features.

 2

		VHL-213 (V1)	p value	VHL-172 (V2)	p value
	Gender
	Male (n=19)	1.48 (0.19-2.84)	0.030	1.07 (0.38-2.73)	0.151
	Female (n=77)	0.81 (0.12-6.48)		0.91(0.25-3.05)
	Age (years)	Corr. Coeff. 0.004	0.971	Corr. Coeff. -0.015	0.885
	Histology
	Classical variant (n=72)	0.93 (0.12-6.48)	0.538	0.87 (0.36-2.73)	0.467
	Other variants (n=24)	0.85 (0.15-4.84)		0.99 (0.25-3.05)
	BRAF
	Wild type (n=39)	1.00 (0.19-3.43)	0.532	0.98 (0.38-2.25)	0.285
	V600E (n=37)	0.93 (0.22-6.48)		0.91 (0.34-3.05)
	Tumor size
	T 1-2 (n=41)	0.94 (0.12-1.96)	0.781	0.94 (0.38-2.25)	0.141
	T 3-4 (n=55)	0.93 (0.15-6.48)		0.91 (0.25-3.05)
	Lymph node metastases
	No (n=57)	0.85 (0.12-6.48)	0.296	0.93 (0.25-2.25)	0.835
	Yes (n=39)	0.99 (0.19-4.84)		0.91 (0.36-3.05)
	TNM stage
	I-II (n=63)	0.8 (0.12-3.43)	0.425	0.86 (0.36-2.73)	0.833
	III-IV (n=33)	0.96 (0.15-6.48)		0.98 (0.25-3.05)
	Recurrences
	No (n=63)	1.04 (0.15-6.48)	0.057	0.94 (0.25-3.05)	0.006
	Yes (n=17)	0.68 (0.22-2.34)		0.63 (0.36-1.6)

  , due to the lack of antibodies able to specifically identify the pVHL 172 .Therefore, we analyzed the mRNA levels of the two VHL splicing variants by means of specific primers. The results obtained indicate that the V1 mRNA variant is the most abundant one in almost all samples analyzed, and the median V1/V2 ratio do not change significantly between tumor and normal matched tissues. The expression of both V1 and V2 is deregulated in the majority of PTC tissues (81.3% and 67.7%, respectively), compared to their normal matched counterparts.Univariate analysis clearly showed the absence of correlation between the VHL mRNAs and patient's age, histological variant, tumor size, lymph node metastases and stage. Reduced V1 mRNA levels significantly associated with female gender, but the meaning of such association remains to be defined. The activating mutation V600E of the BRAF gene is frequently encountered in PTC where it has been associated with a more aggressive phenotype, thus we decided to investigate its effect on the expression of the two VHL mRNA variants. The statistical results

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

immunohistochemistry

Funding

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Author disclosure statement

Authors declare that there is no conflict of interest that could be perceived as prejudicing the impartiality of the research reported.

Highlights

• Reduced expression of the VHL gene has been associated with more aggressive PTCs.

• Here we characterized the expression of the two VHL mRNA variants in PTC tissues.

• We showed that VHL gene expression is deregulated in the majority of PTC tissues.

• Of note is the protective role exerted by VHL transcripts against PTC recurrences.