Isabelle Charpentier

Komlanvi Lampoh

Sensitivity computations in higher order method continuation methods

Keywords: continuation, homotopy, sensitivity, automatic differentiation, Diamant, complex nonlinear eigenvalue problem 2010 MSC: 35B60, 35C20, 65D25, 68W30, 90C31

Sensitivity analysis is a key tool in the study of the relationships between the input parameters of a model and the output solution. Although sensitivity analysis is extensively addressed in the literature, little attention has been brought to the methodological aspects of the sensitivity of nonlinear parametric solutions computed through a continuation technique.

This paper proposes four combinations of sensitivity analysis with continuation and homotopy methods, including sensitivity analysis along solution branches or at a particular point. Theoretical aspects are discussed in the higher order continuation framework Diamant. The sensitivity methods are applied to a thermal ignition problem and some free vibration problems. Remarkable eigenvalue maps are produced for the complex nonlinear eigenvalue problems.

Introduction

Modeling in engineering frequently results in nonlinear parameterized problems, the solutions of which depend on a set of modeling parameters p. These may define material properties, a behavior law, a geometry, initial or boundary conditions. Under-determined nonlinear parametric problems include a particular parameter λ to be varied. Their solutions form branches in bifurcation diagrams [START_REF] Keller | Lectures on numerical methods in bifurcation problems[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. Under analyticity assumptions, the solutions of nonlinear problems may be computed through Taylor-based numerical methods and higher order Automatic Differentiation (AD) [START_REF] Chang | Solving ordinary differential equations using taylor series[END_REF][START_REF] Guillaume | Computation of high order derivatives in optimal shape design[END_REF][START_REF] Berz | COSY INFINITY and its applications in nonlinear dynamics[END_REF][START_REF] Nedialkov | Solving differential-algebraic equations by taylor series (I) : Computing taylor coefficients[END_REF]. Among them, Diamant [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche Diamant[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] is targeted to the solution of nonlinear parametric problems by continuation.

In numerous scientific domains, sensitivity analysis is of prime importance in the understanding of the relationships between the model inputs (modeling parameters, initial or boundary conditions, for instance) and the solution or the branch of solutions. Usages are concerned with uncertainty measurements and optimization, for instance. Sensitivity analysis is thus extensively addressed in the literature [START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF][START_REF] Keulen | Review of options for structural design sensitivity analysis. part 1. linear systems[END_REF]. Nevertheless, to date, little attention has been brought to the methodological study of the sensitivity of nonlinear parametric solutions. In other words, the sensitivity analysis to a perturbation in the modeling parameter is often addressed qualitatively only, by plotting bifurcation diagrams for different values of the parameters.

The paper discusses several quantitative sensitivity techniques adapted to the higher order continuation and homotopy methods proposed by the Diamant framework. Given a scientific problem, the choice for one or another depends on both the continuation objectives (trajectory computations [START_REF] Verschelde | Numerical homotopy algorithms for satellite trajectory control by pole placement[END_REF], path planning [START_REF] Bhattacharya | Search-based path planning with homotopy class constraints[END_REF], complex nonlinear eigenvalue calculation [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF]) and the usage of the sensitivity results. The implementation partially relies on Diamanlab [START_REF] Charpentier | Diamanlab -An interactive Taylor-based continuation tool in MATLAB[END_REF] that takes advantage of the AD generality [START_REF] Griewank | Evaluating Derivatives : Principles and Techniques of Algorithmic Differentiation[END_REF] and provides a graphical user interface for an interactive continuation. The sensitivity drivers are validated on a thermal ignition problem [START_REF] Bazley | The disappearance of criticality in the theory of thermal ignition[END_REF] and the damped beam problem described in the NLEVP benchmark [START_REF] Betcke | NLEVP : A collection of nonlinear eigenvalue problems[END_REF]. A frequency-dependent damping model is then considered to propose remarkable eigenvalue maps for the complex nonlinear eigenvalue problems.

The layout of the paper is as follows. Section 3 discusses the Diamant framework through a brief presentation of both the continuation driver and the homotopy driver. It moreover provides a short introduction to AD and Diamanlab. Theoretical aspects and possible usages of sensitivity computations are described in Section 4 meanwhile numerical results are reported in Section 5. Section 6 provides conclusions.

Setting of the problem

The solutions of a nonlinear parametric problem form branches (u(p), λ(p)) usually represented in a projected bifurcation diagram (Fig. 1) as plots of particular state vector components u(p) against the parameter λ(p). The problem of thermal ignition [START_REF] Bazley | The disappearance of criticality in the theory of thermal ignition[END_REF] is chosen as an illustration. The temperature distribution in a material with exothermic reaction may be written as

∂u ∂t = ∆u + λ exp u 1 + ǫu , (1)
where u = (T -T s )E a /(RT 2 s ) is a normalized temperature variable depending on the temperature T , the surface temperature T s , the universal gas constant R and the Arrhenius activation energy E a . In equation [START_REF] Bazley | The disappearance of criticality in the theory of thermal ignition[END_REF], the parameter ǫ defined as RT s /E determines the normalized temperature behavior. The steady-state problem in the unit square domain,

∆u + λ exp u 1 + ǫu = 0 in Ω = (0, 1) × (0, 1), (2)
u = 0 on the boundary ∂Ω, (3)
is discretized using the five-point stencil finite difference method for the Laplace operator. Eleven equidistributed points are used in each direction. The discrete parametric problem is solved using Taylor-based expansions with a truncation order of 15 and the Diamanlab software (paragraph 3.3.2). The initial guess (u 0 , λ 0 ) is such that u 0 ≡ 0 in Ω and λ 0 = 0. Figure 1 plots the solutions for various values of the parameter ǫ. The case ǫ = 0 corresponds to the Gelfand equation and exhibits a turning point at (u(0.5, 0.5), λ 0 ) ≃ (6.793, 1.385). For 0 < ǫ < ǫ * , the solution branch presents two turning points λ 1 and λ 2 . From a physical point of view, the material ignition takes place when the turning point λ 1 is passed. The solution branch degenerates for ǫ * ≃ 0.24174. This value depends on the geometry of the domain [START_REF] Bazley | The disappearance of criticality in the theory of thermal ignition[END_REF]. Former explanations are those of a brief qualitative sensitivity analysis of the solutions of ( 2)-(3) performed by plotting curves for various values of the parameter ǫ. The quantitative sensitivity analysis of the solution branches requires the computation of the first order derivatives of (u(p), λ(p)) with respect to p in order to propose accurate first order approximations for the solutions.

Computing quantitative sensitivities is an important issue not only for the understanding of the relationships between the nonlinear parametric problem inputs and the output solution branches, but also in classical uses of sensitivity results such as optimization or uncertainty measurement. This paper proposes different options for an accurate calculation of the sensitivities of the solution branches computed through higher-order continuation methods. Sensitivity branches (u d , λ d ) computed using additional Taylor-based expansions are presented in paragraph 4.1. Then, local sensitivities (u p 1 , λ p 1 ) computed at particular points are proposed in paragraph 4.2. Numerical results, including the thermal ignition problem, are reported in Section 5. Prior to that, the next section presents the Diamant framework.

Higher order continuation framework

Diamant [START_REF] Charpentier | Différentiation automatique de la méthode asymptotique numérique typée : l'approche Diamant[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] is the AD-based framework derived from the asymptotic numerical method [START_REF] Thompson | The nonlinear perturbation analysis of discrete structural systems[END_REF][START_REF] Cochelin | Méthode Asymptotique Numérique[END_REF]. Under analyticity assumptions, it automates the solution of the nonlinear parametric problems written in the generic form

R(u(a, p), λ(a, p)) = 0, (4)
a =< u(a, p)u(0, p), ∂u ∂a (0, p) > + (λ(a, p)λ(0, p)) ∂λ ∂a (0, p).

(5)

In system ( 4)-( 5), the residual vector R ∈ C n depends on the unknown state vector u(p) ∈ C n and the scalar parameter λ(p) to be varied. The variable p denotes the modeling parameters. Equation ( 5) measures the pseudo-arc length along the branch of solutions (u(a, p), λ(a, p)) issued from (u(0, p), λ(0, p)). Another path equation may be used to close the under-determined nonlinear residual problem (4). The branch (u(a, p), λ(a, p)) is approximated using the truncated Taylor series [START_REF] Bischof | Adifor -generating derivative codes form fortran programs[END_REF] in a vicinity of a = 0 following

(u(a, p), λ(a, p)) ≃ κ k=0 a k k! ∂ k u ∂a k (0, p), κ k=0 a k k! ∂ k λ ∂a k (0, p) = κ k=0 a k u k , κ k=0 a k λ k , ( 6
)
where κ is the truncation order and u k and λ k denote the Taylor coefficients of u and λ at order k, respectively. The Taylor coefficients of R are denoted by R k and verify

R(u(a, p), λ(a, p)) ≃ κ k=0 a k R k = 0. (7)
For the sake of conciseness, the dependence of u, λ and R with respect to a and p is omitted in the description of the continuation driver (paragraph 3.1) and the homotopy driver (paragraph 3.2). From a theoretical and computer point of view, AD (paragraph 3.3) is a key tool for the evaluation of the Taylor coefficients.

Continuation driver

The analyticity of R implies that

R k = 0, for k = 1 to κ. (8)
The truncated series [START_REF] Bischof | Adifor -generating derivative codes form fortran programs[END_REF] are introduced in the nonlinear problem [START_REF] Bhattacharya | Search-based path planning with homotopy class constraints[END_REF]. Using the Faà di Bruno generalized formula [START_REF] Koutsawa | A generic approach for the solution of nonlinear residual equations. Part I : the Diamant toolbox[END_REF], the Taylor coefficients R k of the compound function R are expressed as

R k = {R 1|u 1 =Id,λ 1 =1 }(u k , λ k ) T + {R k|u k =0,λ k =0 } = 0, for k = 1 to κ, (9)
and comprise the same Jacobian {R 1|u 1 =Id,λ 1 =1 } over the orders. This Jacobian contains the tangent linear contributions resulting from the differentiation of R with respect to u and λ. The remaining terms {R k|u k =0,λ k =0 } contain the nonlinear contributions depending on the Taylor coefficients u l and λ l , for l = 1 to k -1. The resulting κ linear systems,

{R 1|u 1 =Id,λ 1 =1 }(u k , λ k ) T = -{R k|u k =0,λ k =0 }, for k = 1 to κ, ( 10
) a = (u k -u 0 )u 1 + (λ k -λ 0 )λ 1 , for k = 1 to κ. ( 11
)
allow to determine a higher order approximation for the solution (u, λ) of the generic problem (4)- [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF]. The solution of (10)-( 11) is performed in an iterative manner from k = 1 to κ, by alternating higher order derivative computations and linear system solutions. The continuation driver (10)-( 11) is thus a high level mathematical solver. Implementation detail is reported in [START_REF] Charpentier | Diamanlab -An interactive Taylor-based continuation tool in MATLAB[END_REF].

As discussed in previous papers, Diamant is a generic and competitive path following method that may replace the Newton-Raphson scheme in the solution of any equilibrium system of smooth equations. The step size control is based on the a posteriori estimation [START_REF] Charpentier | Exact numerical derivatives of the pair-correlation function of simple liquids using the tangent linear method[END_REF] of the range of validity of the series

a max =        β |u 1 | |u κ |        1/κ-1 , ( 12
)
where β is a small parameter (usually β = 1.10 -6 ). Residuals are measured at the end of the computed branch of solutions to evaluate the interest of a correction. In practice, the higher order predictor Diamant is frequently run without the need for correction steps.

Best practices are concerned with nonlinear mechanical problems involving material or geometrical nonlinearities [START_REF] Cochelin | Méthode Asymptotique Numérique[END_REF][START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF].

Homotopy driver

A homotopy is a continuous deformation that allows for the solution of a nonlinear problem from a known solution using a continuation method [START_REF] Keller | Global homotopies and newton methods[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. Among the best practices relying on homotopy, one may cite trajectory computations [START_REF] Verschelde | Numerical homotopy algorithms for satellite trajectory control by pole placement[END_REF], path planning [START_REF] Bhattacharya | Search-based path planning with homotopy class constraints[END_REF], and complex nonlinear eigenvalue problems [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF].

In the Diamant framework, the homotopy related to the nonlinear residual problem (4) is written [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF] as

R(u(a), λ(a)) = S(u(a), λ(a)) + aT (u(a), λ(a)) = 0, for a ∈ (0, 1). ( 13
)
This homotopy is run from a known solution (u(0), λ(0)) of the problem

S(u(0), λ(0)) = 0, for a = 0, (14)
to the desired solution (u(1), λ(1)) of the residual problem (4), for a = 1. The function T (u, λ) defines as

T (u, λ) = R(u, λ) -S(u, λ), (15)
usually contains the nonlinearities. The use of the Taylor series (6) in equation ( 13) yields generic linear systems

({S 1|u 1 =Id,λ 1 =1 } + a ( j) {T 1|u 1 =Id,λ 1 =1 })(u k , λ k ) T = -{S k|u k =0,λ k =0 } -a ( j) {T k|u k =0,λ k =0 } -T k-1 , ∀ k = 1, .., κ, (16)
that are similar to the linear systems [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF]. In [START_REF] Elkhaldi | A gradient method for viscoelastic behaviour identification of damped sandwich structures[END_REF], the Jacobian related to the j th branch is denoted by J ( j) , that is

J ( j) = ({S 1|u 1 =Id,λ 1 =1 } + a ( j) {T 1|u 1 =Id,λ 1 =1 }). (17)
The linear systems [START_REF] Elkhaldi | A gradient method for viscoelastic behaviour identification of damped sandwich structures[END_REF] are closed with either a path equation or some additional constraint. The choice of S and T and the closure equation is problem-dependent. Other computational aspects may be automated. Initially, this homotopy driver was developed for the solution of the generic complex nonlinear eigenvalue problem [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF] arising in the free vibration modeling of viscoelastic structures [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF]

R(u, λ) = (K + (D 0 + D(λ)) -λM)u = 0, (18)
assuming that the eigensolutions satisfy u(x, t) = u(x)e iωt and ω 2 = λ. In equation ( 18), the stiffness matrix K and the mass matrix M are real symmetric positive definite. The damping term D 0 + D(λ) is a complex symmetric positive semidefinite matrix that may depend on λ in a nonlinear fashion. In [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF], the homotopy is run from one of the eigensolutions of the real problem [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF],

S(u, λ) = (K + D 0 -λM)u, (19)
T (u, λ) = D(λ)u, (20)
to the related eigensolution of the complex problem [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF]. The mass orthonormalization of the eigenvectors

λ k = - u T 0 {S k|u k =0,λ n =0 } + a ( j) {T k|u k =0,λ k =0 } + T k-1 u T 0 ({S 1|u 1 =0,λ 1 =1 } + a ( j) {T 1|u 1 =0,λ 1 =1 }) , (21)
is used to deduce the closure equation.

Automatic differentiation (AD)

Given a simulation code, sensitivity computations and higher order derivatives may be implemented using Automatic Differentiation [START_REF] Griewank | Evaluating Derivatives : Principles and Techniques of Algorithmic Differentiation[END_REF]. AD tools apply the classical differentiation rules such as "the derivative of a product of two functions is the derivative of the first times the second, plus the derivative of the second times the first", and so forth. On the one hand, the source transformation tools [START_REF] Hascoët | TAPENADE 2.1 user's guide[END_REF][START_REF] Bischof | Adifor -generating derivative codes form fortran programs[END_REF] propose to augment the user code with statements encoding the first order derivative computation in tangent linear mode or adjoint mode. The generated code is readable, notably in the tangent linear mode. Some of these tools propose a second order differentiation too. On the other hand, the derivative computations may be hidden in an operator overloading library that includes specific data types as well as methods that encode the classical differentiation rules. The user code is linked to the library at compile time (or at interpretation time) to enable derivative computations at run time. In both cases, the derivative calculation is correct to machine precision.

AD is undoubtedly the more practicable approach for the higher order differentiation, combining generality and ease of use in efficient tools [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF][START_REF] Pryce | AD01, a Fortran 90 code for automatic differentiation[END_REF][START_REF] Charpentier | Fast higher-order derivative tensors with rapsodia[END_REF][START_REF] Berz | COSY INFINITY and its applications in nonlinear dynamics[END_REF][START_REF] Charpentier | Fast higher-order derivative tensors with rapsodia[END_REF]. The higher order AD relies on operator overloading as the vehicle of attaching recurrent derivative computations to the arithmetic operators and intrinsic functions provided by the programming language. A comprehensive overview of the techniques, tools and usage is available on the AD community's website www.autodiff.org. Main applications are sensitivity analysis, gradient-based optimization and nonlinear problem solution.

AD usages in the Diamant continuation framework

Devoted to the generic solution of nonlinear parametric problems, the Diamant drivers rely on the recurrence formulas ( 10) and ( 16), respectively. Each of these sequences of linear systems involves a Jacobian and higher order derivatives. The Jacobian calculation is a well-known issue in AD. It may be achieved using a tangent linear code together with either the canonical basis for small systems, or a sparse Jacobian evaluation for larger ones. Within Diamant, the Jacobian and the higher order terms are computed using an AD operator overloading library. Consequently, Diamant inherits from AD's generality and ease of use.

It is well-known that a particular differentiation is mandatory for an efficient AD of high level mathematical operations such as linear solvers, fast Fourier/wavelet transforms and nonlinear solvers. This also applies to the Diamant continuation drivers.

Diamanlab

Diamanlab [START_REF] Charpentier | Diamanlab -An interactive Taylor-based continuation tool in MATLAB[END_REF] is an object-oriented software in Matlab targeted to the interactive solution of differential algebraic equation systems. It implements the continuation driver [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF]. Solutions branches of the nonlinear problem under study are computed as a set of Taylor-based solutions stored in checkpoints [START_REF] Griewank | Treeverse : An implementation of checkpointing for the reverse or adjoint mode of computational differentiation[END_REF][START_REF] Charpentier | Checkpointing schemes or adjoint codes : Application to the meteorological modelMeso-NH[END_REF]. Each checkpoint contains the Taylor series information (starting point (u(0), λ(0)), series and range of validity) and an oriented tangent vector that indicates the continuation direction. Room may be added to collect complementary information such as sensitivities. During the continuation process, the end point of the last computed branch becomes the next checkpoint unless the user interacts with the plotted bifurcation diagram to set another point and to continue.

From a developer point of view, the core Diamanlab classes are concerned with the handling of the generic nonlinear problem R, the higher-order AD, the checkpoint management and the interactive continuation (graphical user interface). The higher order differentiation methods are gathered in the @T aylor class. The interested reader is referred to [START_REF] Charpentier | Diamanlab -An interactive Taylor-based continuation tool in MATLAB[END_REF] for detail. Diamanlab offers generality and ease of use at the user level too. Any analytical nonlinear satisfying (4) may be solved deriving a user class from the generic system class. For the sake of reproducibility, Diamanlab includes a benchmark of classical nonlinear problems. These may serve as a basis for the implementation of a new user problem.

The sensitivity driver ( 23)-( 24) presented in paragraph 4.1 has been implemented introducing new classes based on the existing Diamanlab classes. The sensitivity driver benefits from the graphical user interface. For now, the homotopy sensitivity driver makes use of the @T aylor class only.

Sensitivity computations and usage

Sensitivity analysis and continuation methods are two of the major pillars in computational sciences and engineering. Designing AD-based numerical methods for accurate quantitative sensitivity computations is then a natural objective. Several options exist. First, sensitivity branches (u d , λ d ) may be approximated using series expanded in the path variable. Second, local sensitivities denoted by (u p 1 , λ p 1 ) may be computed at particular points. These local sensitivities do not depend on the path variable. Both methods are discussed for the continuation and the homotopy drivers, together with their benefits and usages.

Sensitivity branches

The Diamant continuation drivers [START_REF] Charpentier | On higher-order differentiation in nonlinear mechanics[END_REF] and ( 16) are differentiated with respect to the modeling parameter p in order to compute sensitivities as solution branches (u d , λ d ). Such a process involves the Taylor coefficients (u d k , λ d k ),

u d k = ∂u k ∂p = ∂ k+1 u ∂p∂a k (0, p) and λ d k = ∂λ k ∂p = ∂ k+1 λ d ∂p∂a k (0, p), for k = 1, .., κ, (22)
that are mixed derivatives of first order in p and of higher order in a.

Continuation

The tangent linear differentiation of ( 10) with respect to p yields the κ linear systems [START_REF] Lampoh | A generic approach for the solution of nonlinear residual equations. Part III : Sensitivity computations[END_REF] {R

1|u 1 =Id,λ 1 =1 }(u d k , λ d k ) T = -{R d k|u k =0,λ k =0 } -{R d 1|u 1 =Id,λ 1 =1 }(u k , λ k ) T , ∀ k = 1, .., κ, (23)
a d = (u d k -u d 0 )u 1 + (u k -u 0 )u d 1 + (λ d k -λ d 0 )λ 1 + (λ k -λ 0 )λ d 1 , ∀ k = 1, .., κ. (24)
These systems involve the Jacobian exhibited in the continuation driver and may be solved at almost no extra cost. The sensitivity branches (u d , λ d ) are obtained from a differentiation of Eq. ( 6), that is

(u d (a, a d ), λ d (a, a d )) = κ k=0 a k u d k + ka k-1 a d u k , κ k=0 a k λ d k + ka k-1 a d λ k . ( 25
)
From a computer point of view, mixed derivatives [START_REF] Charpentier | Exact numerical derivatives of the pair-correlation function of simple liquids using the tangent linear method[END_REF][START_REF] Charpentier | Sensitivity of solutions computed through the asymptotic numerical method[END_REF][START_REF] Lampoh | A generic approach for the solution of nonlinear residual equations. Part III : Sensitivity computations[END_REF][START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF]] may be computed using a source transformation AD tool onto the code implementing R (the actual user code) to generate a tangent linear code R d . The sensitivity driver implemented using Diamanlab solves equations ( 10)- [START_REF] Charpentier | Diamanlab -An interactive Taylor-based continuation tool in MATLAB[END_REF] as well as ( 23)-( 24). The higher-order AD of R d is performed at interpretation time by defining u d and λ d as @T aylor objects. The higher order tensor algorithm [START_REF] Griewank | Evaluating higher derivative tensors by forward propagation of univariate Taylor series[END_REF], implemented within Rapsodia [START_REF] Charpentier | Fast higher-order derivative tensors with rapsodia[END_REF], is less appropriate for the evaluation of such mixed derivatives because the Diamant drivers alternate derivative computations with linear system solutions.

The sensitivity driver is new to Diamanlab. It is used in paragraph 5.1 to plot branches of sensitivities for the thermal ignition problem. Some sensitivity results obtained for different viscoelastic sandwich beams and plates are reported in [START_REF] Lampoh | A generic approach for the solution of nonlinear residual equations. Part III : Sensitivity computations[END_REF].

Homotopy continuation

The homotopy formulation ( 16) is differentiated with respect to the modeling parameter p following the differentiation process presented in paragraph 4.1.1. The differentiation of ( 16) and (17) yields

J d ( j) = {S d 1|u 1 =Id,λ 1 =1 } + a ( j) {T d 1|u 1 =Id,λ 1 =1 } + a d ( j) {T 1|u 1 =Id,λ 1 =1 }, (26)
J ( j) (u d k , λ d k ) T = -{S d k|u k =0,λ k =0 } -a ( j) {T d k|u k =0,λ k =0 } -a d ( j) {T k|u k =0,λ k =0 } -T d k-1 -J d ( j) (u k , λ k ) T , ∀ k = 1, .., κ. (27)
The closure equation is to be differentiated too. The sensitivity of solution branches is of special interest when the trajectory computation is the main objective of the homotopy [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF][START_REF] Verschelde | Numerical homotopy algorithms for satellite trajectory control by pole placement[END_REF], including the evolution of the complex eigenvalues with respect to a particular modeling parameter as proposed in paragraph 5.2.2.

Local sensitivities

Sensitivities along branches may have limited interest when the goal of the homotopy is the solution at the end point only. Likewise, the sensitivity of solution branches computed and plotted using an interactive continuation software may be questionable when the user operates a "jump" from one point of the bifurcation diagram to another, for instance. The sensitivity drivers should thus allow for a local sensitivity analysis, at either a user-defined checkpoint or at the end point of the homotopy method. In both cases, the differentiation stages in a and p are decoupled and there are no more mixed derivatives to compute.

Continuation

Let (u(p), λ(p)) be a known solution point. An approximated solution of (4) may be computed for a small perturbation of p in the direction of the unit vector p 1 as

(u(p + αp 1 ), λ(p + αp 1 )) = (u(p) + αu p 1 (p)), λ(p) + λ p 1 (p)), ( 28
)
where α is a small real parameter and the sensitivity (u p 1 (p), λ p 1 (p)) is the first order derivative of (u(p), λ(p)) computed with respect to p. From a numerical point of view, one may use either the tangent linear code R d , or the code R overloaded for a differentiation with respect to p.

At a given point, local sensitivities do not depend on the path variable a and the path equation ( 5) cannot be used to close the system. As in classical parameterized problems, the local sensitivity of u may be computed for a given λ. For the sake of completeness, it is important to mention that higher-order local approximations may be obtained following the asymptotic analysis proposed in [START_REF] Guillaume | Computation of high order derivatives in optimal shape design[END_REF].

Homotopy

The differentiation of the homotopy and the propagation of the derivatives along the path are not mandatory when one is interested in the sensitivity of the end point solution. The particular case of the nonlinear eigenvalue problem with an orthonormality closure constraint ( 21) is discussed in paragraph 5.2. This local computation requires one extra linear system solution at the end point, but it leaves the Taylor series capabilities available for an optimization method [START_REF] Elkhaldi | A gradient method for viscoelastic behaviour identification of damped sandwich structures[END_REF][START_REF] Guillaume | Computation of high order derivatives in optimal shape design[END_REF][START_REF] Guillaume | Nonlinear eigenproblems[END_REF] based on a higher order differentiation of the eigenvariables u and λ with respect to the modeling parameter p. From a numerical point of view, the paper focuses on sensitivity results only.

Numerical results and discussion

The sensitivity drivers are evaluated on a thermal ignition problem [START_REF] Bazley | The disappearance of criticality in the theory of thermal ignition[END_REF] for the continuation one, and on the damped beam problem [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] described in the NLEVP benchmark [START_REF] Betcke | NLEVP : A collection of nonlinear eigenvalue problems[END_REF] for the homotopy driver. A frequency-dependent damping model [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF][START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF] is considered to emphasize the abilities of the Diamant driver and to discuss approximated eigenvalue maps.

Thermal ignition problem

Branch and checkpoint sensitivity results are proposed for the steady-state problem presented in Section 2. The code related to the discrete thermal ignition problem is denoted by R.

∆u d + λ d + λ u d -u 2 ǫ d (1 + ǫu) 2 exp u 1 + ǫu = 0 in Ω = (0, 1) × (0, 1), (29)
u d = 0 on the boundary ∂Ω.

Equations ( 29)-( 30) are discretized with the five-point stencil finite difference scheme used for the equations ( 2)-( 3) and implemented in the tangent linear code denoted by R d . This contains the statements of R as well. In Matlab, the code R d is overloaded at interpretation time for the higher AD. The sensitivity driver solves the equations ( 10)-( 11) and ( 23)-( 24) to compute the Taylor coefficients and the branches of sensitivities following [START_REF] Hershberger | Computing minimum length paths of a given homotopy class[END_REF]. For the sake of reproducibility, the sensitivity driver and the user class related to the thermal ignition problem may be obtained from the corresponding author. The correctness of the sensitivity process is checked with a classical Taylor test. This compares the sensitivity result (u d , λ d ) obtained with the sensitivity driver at point ǫ in the direction ǫ d to first order finite difference approximations computed using the continuation driver as follows

r α (u) = |u(ǫ + αǫ d ) -u(ǫ)| α|u d (ǫ, ǫ d )| and r α = (λ) |λ(ǫ + αǫ d ) -u(ǫ)| α|λ d (ǫ, ǫ d )| , for small α ∈ I R. ( 31
)
α r α (u) r α (λ) 1.E-3 1.13957691646213 0.99771091315504 1.E-4 1.00808650725132 0.99988326634830 1.E-5 1.00075453289025 0.99998938980647 1.E-6 1.00007491725486 0.99999894943411 1.E-7 1.00000748384382 0.99999989502055 1.E-8 1.00000072256463 0.99999997858741 1.E-9 0.99999752919078 1.00000005788443 1.E-10 0.99999122129183 0.99999883793030 1.E-11 1.00003853053399 0.99999639802205 Table 1: Taylor test on the sensitivity of (u(0.5, 0.5), λ) with respect to (ǫ, ǫ d ) = (0.2, 0.1) at the end point of the 5 th branch. Ratios computed at the end point of the fifth branch using (ǫ, ǫ d ) = (0.2, 0.1) are reported in Table 1. Theoretically, these ratio should tend linearly toward 1 as α tends to 0 since small α are necessary to minimize the numerical errors coming from the finite difference approximation. However, the subtraction of too close floating-point numbers yields a cancellation error that dominates the truncation error coming from the finite difference computation. Taylor tests presenting such a behavior indicate that the implemented sensitivity computation is correct.

Bifurcation diagrams and their first order approximations are displayed in Fig. 2. Solid lines plot solutions (u, λ) computed for various values of the modeling parameter ǫ (blue for ǫ = 0.2, green for ǫ = ǫ * and red for ǫ = 0.3) using the continuation driver. The solution (u d , λ d ) is calculated using the sensitivity driver with the parameter ǫ = 0.2 and the perturbation ǫ d = 0.1. The first order approximated solutions are plotted using dashed lines with same colors. The agreement between the thermal ignition problem solutions and their approximations using the sensitivity computations is good. In particular, the turning point positions are well-reproduced for this highly nonlinear problem (Figs. 2 and3.a).

In Fig. 3.a, the sensitivity (u d , λ d ) evaluated along the path for (ǫ, ǫ d ) = (0.2, 0.02) is represented as blue dotted lines linking the solution computed at ǫ = 0.2 (blue curve) to the first order approximation (u + u d , λ + λ d ) (magenta dashed curve). The solution computed at ǫ = 0.22 is plotted using a magenta solid line. Checkpoints are indicated using circles with same colors. This picture shows that the proposed continuation driver requires very few branches computations, about twenty, each one involving a Jacobian computation. The first order approximation of the solution at ǫ = 0.22 is very good. One notices that the checkpoint locations slightly differ due to the adaptive computation of the validity range a max . This is not important in the proposed numerical method. Nevertheless, this may be a significant issue in the sensitivity computations performed by means of finite difference approximations.

Local sensitivity branches, denoted by (u+u d , λ+λ d ) too, are evidenced in Fig. 3.b. These are computed initializing (u d 0 , λ d 0 ) to (0,0) at checkpoints to cancel the linear accumulation in the first order derivatives.

Nonlinear eigenvalue problems

Initially, the Diamant homotopy driver is designed for the solution of the complex frequency-dependent eigenvalue problems arising in the free vibration modeling [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF]. For the sake of reproducibility, sensitivity computations are illustrated on the damped_beam problem [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] belonging to the collection of nonlinear eigenvalue problems NLEVP [START_REF] Betcke | NLEVP : A collection of nonlinear eigenvalue problems[END_REF].

NLEVP's damped beam example

The damped_beam problem [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] is a quadratic eigenvalue problem (QEP) [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] arising in the free vibration analysis of the simply supported beam, Fig. 4, with a damper located in its middle. The damping coefficient d c is constant. Assuming that the eigensolutions satisfy u(x, t) = u(x)e ωt , the resulting discrete damped_beam QEP is written as

Q(ω)u = (ω 2 M + ωD + K)u = 0, ( 32
)
where the stiffness matrix K and the mass matrix M are symmetric positive definite, and the viscosity matrix D is symmetric positive semidefinite. The damped_beam modeling parameters are reported in Tab. 2. In NLEVP, the quadratic polynomial Q(ω) is turned into a linear one,

L(ω) = ωX + Y = 0, ( 33
) Length L = 1 m Young modulus E = 7 × 10 10 Pa Width l = 5 × 10 -2 m Density ρ = 0.674 Kg.m -3 Height H = 5 × 10 -3 m Damping coefficient d c = 5
Table 2: Geometric and material properties of the damped beam.

-16 - the dimension of which is twice those of Q(ω). The damped_beam raw eigenvalues (Fig. 5.a, no scaling is applied) are obtained concurrently using a QZ algorithm following where the cell array coeffs={K,D,M} contains the problem matrices.

The authors make use of this example to draw attention on important linearization aspects such as scaling, condition numbers and backward errors. Accurate results (Fig. 5.b) may be obtained either by scaling the QEP before the solution using a linearization [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF], or by using the homotopy driver as proposed in the next paragraph.

Diamant & homotopy

Substituting ω by i √ λ into the QEP (32) yields the residual problem satisfying (18)

(K + i √ λD -λM)u = 0, (34)
where K, M and D are the damped_beam matrices. The problem ( 34) is turned into the generic complex nonlinear eigenvalue problem [START_REF] Griewank | ADOL-C, a package for the automatic differentiation of algorithms written in C/C++[END_REF] choosing D 0 = 0 and D(λ) = i √ λD. The complex eigensolutions of (34) are determined from the real ones in a sequential manner. A unique branch computation is necessary per eigenvalue. The eigenvalues with positive imaginary parts are computed with Diamant and plotted in Fig. 5.b. The others are reproduced by symmetry. The computed spectrum agrees with those obtained in [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] for the best linearization.

The Matlab version of the Diamant homotopy driver is obviously more time consuming than NLEVP in the solution of [START_REF] Nayak | Dynamic stability of magnetorheological elastomer based adaptive sandwich beam with conductive skins using FEM and the harmonic balance method[END_REF]. Nevertheless, the homotopy driver has a number of precious methodological assets. First, the automation and generality of Diamant and AD allows for the solution of free vibration problems with nontrivial damping laws [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF] such as a generalized Maxwell model or a Kelvin-Voigt model. Second, the Taylor-based approximations are obtained with an a posteriori error measured and controlled along the path. The linearization phase discussed in [START_REF] Higham | Scaling, sensitivity and stability in the numerical solution of quadratic eigenvalue problems[END_REF] is not necessary. Third, Taylor-based approximations provide access to the continuous variation of each of the eigenvalues along its homotopy path. For instance, Fig. 6 accounts for the evolution of the spectrum with respect to the damping coefficient ad c for a ∈ {0, 0.2, 0.6, 1}. In other words, the homotopy intrinsically provides a sensitivity analysis with respect to the damping function since it is scaled by the path parameter a.

For illustration purposes, the damping is now achieved considering the frequency-dependent rheological model determined for the viscoelastic material 3M ISD112 at a temperature of 27 • [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF],

G(ω) = G 0 1 + 3 j=1 ∆ j ω ω -iΩ j , (35)
where G 0 is the shear modulus of the delayed elasticity and ω = √ λ. Parameters (∆ j , Ω j ) are reported in Tab. 3. The damping is still applied at the middle of the beam using the matrix D = D/d c . As in [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF], the real constant part D 0 and frequency-dependent part D(ω) of the damping function (including the height of H = 5 × 10 -3 m),

D 0 + D(ω) = HG(ω) D, (36)
are split into the homotopy functions S and T as follows S(u, λ) = (K + HG 0 D -λM)u, and T (u, λ) = HG 0 3 j=1

∆ j ω ω -iΩ j Du. (37)
Nonlinear mechanics results obtained at homotopy end points are discussed in [START_REF] Bilasse | A generic approach for the solution of nonlinear residual equations. Part II : Homotopy and complex nonlinear eigenvalue method[END_REF] for viscoelastic sandwich beams and plates. The splitting [START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF] may also be used,

S(u, λ) = (K -λM)u, and T (u, λ) = HG 0 1 + 3 j=1 ∆ j ω ω -iΩ j Du. ( 38
)
The choice between these two splittings may be of particular interest when an electro-rheological law [START_REF] Yeh | Finite element dynamic analysis of orthotropic sandwich plates with an electrorheological fluid core layer[END_REF] or a magneto-rheological law [START_REF] Nayak | Dynamic stability of magnetorheological elastomer based adaptive sandwich beam with conductive skins using FEM and the harmonic balance method[END_REF] is used. The difference between these models and viscoelastic ones is in the parameter E * that represents an applied electrical or magnetic field to be varied. For instance, the shear modulus used in [START_REF] Yeh | Finite element dynamic analysis of orthotropic sandwich plates with an electrorheological fluid core layer[END_REF] reads as

G(E * ) = m g E * + b g , (39)
Figure 7: Eigensolutions for various scaling of the damping model [START_REF] Thompson | The nonlinear perturbation analysis of discrete structural systems[END_REF]. (a) : Splitting [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF]. (b) : Splitting [START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF] where the parameters m g and b g are obtained from experiments. It does not depends on ω. A homotopy based on the splitting [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF], with T = m g E * , enables for the observation of the evolution of the eigensolutions with respect to the applied electrical field. Work is in progress to consider such a law in an actual mechanical framework. Meanwhile, Fig. 7 plots the eigenvalue evolution for the beam damped with the 3M ISD112 function considering the splittings [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF] and [START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF], respectively. A unique branch computation is necessary per eigenvalue, except for the first eigenvalue computed with the splitting [START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF].

The half spectrums (positive imaginary parts) corresponding to damping function (D 0 +D(ω)) for various a ∈ (0, 1) are plotted in Figs. [START_REF] Chang | Solving ordinary differential equations using taylor series[END_REF].a and 7.b for the splitting [START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF] and [START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF], respectively. The choice for one or another splitting influences the construction of the real eigenvalue problem and determines the "slope" between the real eigenvalue (starting point of the homotopy) and the related complex eigenvalue (end point of the homotopy). Intermediate results are different along the path since the path parameter a acts as a scaling on T only. As expected, the first eigenvalue is the more sensitive to the splitting choice and both splittings produce equivalent eigensolutions at the end points. The first five complex eigenvalues computed at the end point of the homotopy are reported in Tab. 5 for the two splittings.

Sensitivities

The few statements of the damped_beam code are differentiated (manually) with respect to modeling parameters to get the tangent linear matrices K d , D d and M d . Taylor tests, see Tab. 4, are verified in a systematic manner. Figure 8 presents sensitivity results for a perturbation of the density coefficient ρ and the two splittings. As expected, the sensitivity on ω is negative because a denser beam vibrates slower and has lower frequencies. Computations are performed using ρ = 0.674Kg.m -3 and ρ d = 0.1Kg.m -3 . In Figure 9.a, the eigenvalues (color circles, blue for a = 1, magenta for a = 0.8, green for a = 0.6, cyan for a = 0.4, red for a = 0.2 and black for a = 0) and their related sensitivity are added to propose a first order approximation (black crosses) of the first nine eigenvalues computed with ρ = 0.774Kg.m -3 (with squares and same colors). For the two splittings, the good agreement between squares and crosses is confirmed by relative errors (about 0.8%) presented in Tab. 5 on the first order approximations using the - sensitivity ω d computed along the path.

α r α (λ) r α (λ) (nlevp, d c = 5i) (ISD112, H = 5.10 -3 m)
Sensitivities computed along the path with respect to both the homotopy parameter and one of the modeling parameter, allow for the construction of approximated eigenvalue maps. Similar maps may be constructed for the other modeling parameters. Consequently, approximated eigenvalues may be evaluated in any non trivial direction of perturbation αp 1 (α ∈ I R small, |p 1 | = 1) from the homotopy checkpoints.

At the end point of the homotopy, see paragraph 4.2.2, the local sensitivity ω p 1 of a complex eigenvalue ω may be obtained from a differentiation with respect to p of both the mass orthonormalization formula and the generic equation ( 4). The tangent linear equation ( 40)

Mu p 1 = -0.5M 1 u, (40)
is derived the orthonormalization formula and solved from the known complex eigensolution (u, λ) to compute the sensitivity u It may be computed as follows

λ p 1 = - u T {R 1|λ p 1 =0,u p 1 ,p 1 } u T {R 1|λ p 1 =1,u p 1 =0,p 1 =0 } , (42)
to deduce

ω p 1 = 1 2 λ p 1 √ λ . (43)
In equation ( 42), the tangent linear term {R 1|λ p 1 =1,u p 1 =0,p 1 =0 } is the first order derivative of R computed with respect to λ only. The complementary term {R 1|λ p 1 =0,u p 1 ,p 1 } contains the first order derivatives with respect to p and u evaluated using the direction of perturbation p 1 and the sensitivity u p 1 . Thinking in further higher order approximations, computations are performed using the code R overloaded for a differentiation with respect to p. The tangent linear code R d may also be used. Both methods are implemented and quasi-identical results are obtained. Some relative errors measured between the perturbed eigensolutions ω(0.774) and their the first order approximation obtained with (41)-(42) are reported in Tab. 5. The different approaches proposed to compute the eigenvalue sensitivities yield very similar results at the end points.

Conclusions

On the one hand, continuation and homotopy are classical tools for the solution of nonlinear parametric problems arising in engineering applications. On the other hand, sensitivity analysis are precious tools for model improvement, conception or uncertainty measurement. Proposing sensitivity drivers for the continuation methods is a natural issue. This paper discusses several options for the sensitivity analysis with respect to the modeling parameters of the solution branches produced by continuation, including homotopy, together with particular applications and usages. Theoretical developments and implementation are realized in the Diamant framework that combines higher order Taylor approximations to AD, inheriting from the generality of Diamant and AD. Numerical results are presented for a thermal ignition problem and the damped beam problem issued from the NLEVP collection. In both cases, accurate sensitivities are computed.

A particular attention is brought to the complex nonlinear eigenvalue problems arising in the free vibration modeling of damped structures. A frequency-dependent rheological model is considered to emphasize the interest of sensitivity computations performed along the homotopy path. More precisely, the sensitivity driver we propose allows for the construction of approximated eigenvalue maps from both a very few Taylor series expanded in the homotopy path parameter and a sensitivity analysis performed concurrently with respect to another modeling parameter. Future work will be concerned with a higher order conception method taking advantage of eigenvalue maps.

Références

Figure 1 :

 1 Figure 1: Bifurcation diagrams, u(0.5, 0.5) against λ, for the thermal ignition problem.

  The sensitivity (u d , λ d ) = ∂u ∂ǫ , ∂λ ∂ǫ of the solution (u, λ) of the thermal ignition model (2)-(3) satisfies

Figure 2 :

 2 Figure 2: Bifurcation diagrams, u(0.5, 0.5) against λ, and their first order approximations for the thermal ignition problem. The sensitivity is computed at ǫ = 0.2. Color online.

Figure 3 :Figure 4 :

 34 Figure 3: Sensitivities of the solution computed at ǫ = 0.2. (a) Along the path, ǫ d = 0.1ǫ. (b) : Local branches at checkpoints. Color online.L

Figure 5 :

 5 Figure 5: Eigensolutions of the NLEVP damped beam problem. (a) : NLEVP raw results. (b) : Homotopy results.

[

  coeffs, fun] = damped_beam(200); [X, e] = polyeig(coeffs{:});

Figure 6 :

 6 Figure 6: Eigensolutions for various values of the constant damping coefficient d c .

Figure 8 :

 8 Figure 8: Spectrum sensitivity, with respect to ρ, for various values of the damping coefficient. (Color online).

Figure 9 :

 9 Figure9: Eigenvalue maps and first order approximations for a = 1 (blue), a = 0.8 (magenta), a = 0.6 (green), a = 0.4 (cyan), a = 0.2 (red), a = 0 (black). (a) : Homotopy (37) ; (b) : Homotopy[START_REF] Utke | Sensitivity analysis for mixed-language numerical models[END_REF]

p 1 of 1 + {R 1|λ p 1 =0,u p 1 ,p 1 }

 11111 the eigenvector u. The sensitivity λ p 1 of the eigenvalue verifies {R 1|λ p 1 =1,u p 1 =0,p 1 =0 }λ p = 0.(41)

Table 3 :

 3 Parameters for the 3M ISD112 model at a temperature of 27 • C[START_REF] Trindade | Modeling of frequeney-dependent viscoelastic materials for active-passive vibration damping[END_REF].

	j	∆ j Ω j (rad/s)
	1	0.746	468.7
	2	3.265	4742.4
	3 43.284	71532.5

Table 4 :

 4 Taylor test on the sensitivity of the first complex eigenvalue λ (at homotopy end point) with respect to a perturbation of the density (ρ, ρ d ) = (0.2, 0.1).

	1.E-2 0.9889713342713	0.98895850930239
	1.E-3 0.9989467373101	0.99887758359406
	1.E-4 1.0011346796444	1.00115275243155
	1.E-5 1.0145292086048	1.01473807301412
	1.E-6 1.0772568674870	1.07828544963728
	1.E-7 2.6064844849179	2.61097302946291

Table 5 :

 5 Eigenvalues ω and sensitivities ω d to a perturbation of the density computed at the end point of the homotopy path. Relative errors for the first five eigenvalues and the two splittings.

		Eigenvalue (ǫ = 0.674)	Sensitivity	Eigenvalue (ǫ = 0.774)	Path comp.	End-point
	N o	ω(0.674)	ω d (0.674, .1)	ω(0.774)	Rel. err.	Rel. err.
	Splitting (37)
	1	1.130e+002 + 1.018e+001i -8.539e+000 + -1.496e+000i 1.053e+002 + 8.872e+000i 8.354e-003 8.241e-003
	2	2.904e+002 + 6.424e-018i	-2.154e+001 + -8.358e-019i	2.709e+002 + 4.641e-016i	7.878e-003 7.878e-003
	3	6.622e+002 + 7.035e+000i	-4.932e+001 + -8.365e-001i	6.178e+002 + 6.296e+000i 7.923e-003 7.936e-003
	4	1.161e+003 + 5.077e-021i	-8.616e+001 + -6.396e-022i	1.084e+003 + 2.142e-019i	7.878e-003 7.878e-003
	5	1.819e+003 + 4.916e+000i	-1.351e+002 + -6.365e-001i	1.697e+003 + 4.357e+000i 7.891e-003 7.893e-003
	Splitting (38)
	1	1.130e+002 + 1.027e+001i -8.530e+000 + -1.511e+000i 1.053e+002 + 8.951e+000i 8.341e-003 8.287e-003
	2	2.904e+002 + 3.445e-016i	-2.154e+001 + -4.552e-017i	2.709e+002 + 2.218e-016i	7.878e-003 7.878e-003
	3	6.622e+002 + 7.035e+000i	-4.931e+001 + -8.379e-001i	6.178e+002 + 6.296e+000i 7.921e-003 7.936e-003
	4	1.161e+003 + 1.444e-018i	-8.616e+001 + -1.821e-019i	1.084e+003 + 2.175e-021i	7.878e-003 7.878e-003
	5	1.819e+003 + 4.916e+000i	-1.351e+002 + -6.367e-001i	1.697e+003 + 4.357e+000i 7.891e-003 7.893e-003


Alice Talpin

Félicie Costantino

Nelly Bonilla

Ariane Leboime

Franck Letourneur

Sébastien Jacques

Florent Dumont

Sonia Amraoui

Charles-Antoine Dutertre

Henri-Jean Garchon

Maxime Breban

Gilles Chiocchia
email: gilles.chiocchia@inserm.fr

Monocyte-derived dendritic cells from HLA-B27+ axial spondyloarthritis (

Introduction

Spondyloarthritis (SpA) is a chronic inflammatory rheumatic disorder, with a prevalence of around 0.42% in Caucasian populations [START_REF] Costantino | Prevalence of spondyloarthritis in reference to HLA-B27 in the French population: results of the GAZEL cohort[END_REF]. Joint inflammation is responsible for pain and stiffness but long-term outcome is mainly determined by new bone formation, which can lead to complete ankylosis. Effective treatments such as TNF blockers can temporarily suppress inflammation, but none has yet proven to affect long-term disease outcome. Therefore, there is a need to better understand pathological mechanisms controlling both initiation and progression of SpA.

Genome-wide gene expression analysis is a powerful approach to identify molecular mechanisms responsible for a disease. Microarray studies have already been conducted in SpA, most of them focusing on the ankylosing spondylitis (AS) subtype [START_REF] Smith | Gene expression analysis of macrophages derived from ankylosing spondylitis patients reveals interferon-gamma dysregulation[END_REF][START_REF] Gu | Identification of RGS1 as a candidate biomarker for undifferentiated spondylarthritis by genome-wide expression profiling and real-time polymerase chain reaction[END_REF][START_REF] Duan | Gene expression profiling reveals a downregulation in immune-associated genes in patients with AS[END_REF][START_REF] Assassi | Whole-blood gene expression profiling in ankylosing spondylitis shows upregulation of toll-like receptor 4 and 5[END_REF][START_REF] Pimentel-Santos | Whole blood transcriptional profiling in ankylosing spondylitis identifies novel candidate genes that might contribute to the inflammatory and tissue-destructive disease aspects[END_REF][START_REF] Sharma | Insights in to the pathogenesis of axial spondyloarthropathy based on gene expression profiles[END_REF]. Each of those studies identified several genes differentially expressed between patients and controls but there was very little overlap between their results, maybe due to the variety of protocols.

Several studies support an important role of dendritic cells (DCs) in the pathogenesis of SpA. Indeed, aberrant functions of DCs have been demonstrated in an HLA-B27/human β2-microglobulin transgenic rat model of SpA [START_REF] Breban | Animal models of HLA-B27associated diseases[END_REF]. First, DCs from these rats have a decreased capacity to stimulate primary allogeneic or syngeneic T cell responses [START_REF] Hacquard-Bouder | Defective costimulatory function is a striking feature of antigen-presenting cells in an HLA-B27-transgenic rat model of spondylarthropathy[END_REF]. Furthermore, the proportion of conjugates formed between HLA-B27 DCs and naïve CD4 + T cells is reduced. Moreover, mature HLA-B27 molecules expressed by DCs appear to impair the formation of an antigen-independent immunologic synapse with naïve CD4 + T cells by interfering with the engagement of co-stimulatory molecules [START_REF] Hacquard-Bouder | Alteration of antigen-independent immunologic synapse formation between dendritic cells from HLA-B27-transgenic rats and CD4+ T cells: selective impairment of costimulatory molecule engagement by mature HLA-B27[END_REF]. Finally, B27 transgenic rat DCs were shown to favor Th17 expansion [START_REF] Glatigny | Proinflammatory Th17 cells are expanded and induced by dendritic cells in spondylarthritis-prone HLA-B27-transgenic rats[END_REF][START_REF] Utriainen | Expression of HLA-B27 causes loss of migratory dendritic cells in a rat model of spondylarthritis[END_REF] and to alter regulatory T cell function, resulting in decreased IL-10 and enhanced IL-17 production [START_REF] Araujo | Increased production of interleukin-17 over interleukin-10 by regulatory T cells implicates ICOS molecule in experimental spondyloarthritis[END_REF].

On this basis, we hypothesized that DCs also play an important role in human disease, as proposed in the HLA-B27 transgenic rat model. The aim of the current study was to compare functional capacity and the gene expression profile of monocyte-derived DCs (MD-DCs) in patients with HLA-B27 + axial SpA and healthy controls.

Methods

Patients and controls

We studied three different groups of patients and controls: 10 patients and 5 controls for monocyte subset characterization by flow cytometry, 19 patients and 24 controls for mixed lymphocyte reaction and 9 patients and 10 controls for the MD-DC transcriptomic study (6 patients and 4 controls were investigated for the two latter). All the patients were HLA-B27 + and fulfilled the Assessment of SpondyloArthritis International Society classification criteria for axial SpA [START_REF] Rudwaleit | The development of Assessment of SpondyloArthritis international Society classification criteria for axial spondyloarthritis (part II): validation and final selection[END_REF]. Controls were healthy blood donors.

Patients' and healthy controls' characteristics are summarized in the Additional files (see Additional file 1: Table S1 and Additional file 2: Table S2, respectively). All participants in the study gave written informed consent and the study was approved by local ethics committee of Ile-de-France XI (Saint-Germaine-en-Laye France).

Cell isolation, culture and stimulation

Peripheral blood mononuclear cells (PBMCs) were isolated from 50 mL of blood by gradient separation on Ficoll density gradient centrifugation (STEMCELL Technologies, Grenoble, France). Monocytes used to generate MD-DCs were purified by magnetic cell sorting using anti-CD14 monoclonal antibody (mAb)-coated beads (BD IMag, Le Pont de Claix, France). Sorted monocytes were morphologically homogeneous with 99% of CD14 + cells, as determined by flow cytometry.

Monocytes were further cultured for 6 days in 24-well plates (400,000 cells/500 μL) in Roswell Park Memorial Institute (RPMI) 1640 medium supplemented with 10% heat-inactivated fetal calf serum, 100 U/mL penicillin, 100 μg/mL streptomycin, 500 U/mL recombinant human granulocyte-macrophage colony stimulating factor (rhGM-CSF) and 500 U/mL rhIL-4 (AbCys, Paris, France). Then, the MD-DCs were stimulated or not with lipopolysaccharide (LPS) from Escherichia coli (LPS, Sigma-Aldrich, St Louis, MO, USA) at a concentration of 100 ng/mL for the last 6 or 24 hours of culture (further referred to as time points H0, H6 and H24).

CD4 + T cells were purified from PBMCs from two unrelated healthy donors by magnetic cell sorting using anti-CD4 monoclonal antibody (mAb)-coated beads (BD IMag), and stored frozen until used for mixed lymphocyte reaction (MLR).

Flow cytometry

To characterize monocyte subsets, freshly purified PBMCs were analyzed by six-color flow cytometry on FACS LSRII apparatus. The gating strategy was based on a previous report [START_REF] Dutertre | Pivotal role of M-DC8 + monocytes from viremic HIV-infected patients in TNFα overproduction in response to microbial products[END_REF]. Monocytes were subdivided into three major subsets: classical CD14 ++ CD16 -, intermediate CD14 ++ CD16 + and non-classical CD14 + CD16 ++ monocytes. The following anti-human mAbs were used: CD45-Amcyan (BD Biosciences), HLA-DR-PerCP (BD Biosciences), CD19-ECD (Beckman Coulter), CD14-QDot655 (Invitrogen), CD16-APC-H7 (Beckman Coulter, Villepinte, France). The Live/ Dead blue Dye (Invitrogen) was used to exclude dead cells.

Samples of the purified monocytes used to generate MD-DCs and of the resulting MD-DCs were routinely stained with the following anti-human mAbs: CD14-FITC, CD11c-APC, CD40-PE, HLA-I-FITC, HLA-DR-PerCP, CD80-PE, CD83-APC and CD86-FITC (all from BD Bioscience) and analyzed by flow cytometry on FACS canto II apparatus (BD Biosciences).

Mixed lymphocyte reaction (MLR)

Purified allogeneic CD4 + T cells (10 5 cells per well) from healthy donors were cultured with unstimulated (H0) or LPS-stimulated (H6, H24) MD-DCs (10 4 cells per well), in 96-well flat-bottomed culture dishes in a final volume of 200 μL. Proliferation of T cells was assayed by measuring incorporation of 3 H-deoxythymidine added (0.5 μCi per well) after 6 days of culture, using a Microbeta scintillation counter (Wallac, Turku, Finland). Data are expressed as the mean counts per minute (CPM) in triplicate wells. An MLR index (ratio of CPM of MLR on CPM of CD4 + T cells only) was used to represent CD4 + T cell proliferation. Two stored CD4 + T cell batches from different healthy donors were sequentially used for MLR in two sets of experiments, each including equivalent numbers of patient and control MD-DC samples. As there was no statistically significant difference in the results between both sets of experiments, we pooled them. The Wilcoxon test was used to compare MLR indices between patients and controls at each stimulation time point.

Transcriptomic study RNA isolation

MD-DCs were disrupted and homogenized using RLT buffer (Qiagen, Valencia, CA, USA). Total RNA was isolated using RNeasy Mini Kit (Qiagen). RNA quantity and quality were assessed using Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA, USA). Only samples with an RNA integrity number (RIN) above 8 were further processed.

Microarray hybridization

RNA was reverse-transcribed, converted to biotinylated complementary RNA using standard Affymetrix protocol (Affymetrix, Santa Clara, CA, USA) and hybridized to the Affymetrix GeneChip Human Gene 1.0 ST Array by the genomic platform of the Cochin Institute.

Differential gene expression validation by qRT-PCR

For validation, the relative gene expression levels of candidate genes identified through the foregoing microarray study were further quantified using qRT-PCR. Briefly, RNA treated with DNase I (Invitrogen) was reversetranscribed using SuperscriptII (Invitrogen) and then quantified using the SYBR green PCR Master Mix (Applied Biosystems) and the 7300 Real-Time PCR System (Applied Biosystems). Primers were purchased from Eurofins MWG (nucleotide sequences of the PCR primers are available in Additional file 3: Table S3). The experiment design included three technical replicates.

Statistical analysis

Raw Affymetrix data (.cel files) from 57 arrays (corresponding to H0, H6 and H24 time points of stimulation by LPS for 19 subjects) were transformed by the Robust Multichip Analysis (RMA) method using Bioconductor in R software (library Affy) [START_REF] Gautier | affy-analysis of Affymetrix GeneChip data at the probe level[END_REF]. This transformation included background correction, normalization and summarization of expression values using Brainarray version 15 custom chip definition (cdf) files generated with the Ensembl annotation set [START_REF] Dai | Evolving gene/transcript definitions significantly alter the interpretation of GeneChip data[END_REF]. Before analysis, genes were filtered on annotation and expression level mean with a cut off at 6.11 (corresponding to the mean of expression level of both antigenomic and intronic probes of housekeeping genes). Gene expression levels were then fitted to a bivariate linear model including disease status and time point after stimulation, using the LInear Models for Microarray data (LIMMA) package of Bioconductor in R environment to determine differentially expressed genes [START_REF] Smyth | Linear models and empirical bayes methods for assessing differential expression in microarray experiments[END_REF]. The threshold for differential gene expression between patients and controls was set to a global fold change above or below 1.5 with a nominal P-value below 0.01.

For qRT-PCR validation, gene expression data were computed with the ΔΔCq method, using three housekeeping genes: RPL30, β-ACT and GAPDH. Then, twoway analysis of variance (ANOVA) was used to test differential gene expression among samples from SpA and controls without (baseline; no treatment) and with LPS treatment for 6 and 24 hours. In case of statistical significance (P <0.05), post hoc t-tests with the Bonferroni correction allowed us to perform pairwise comparisons.

Paired genes co-expression study and molecular pathway analysis

Pairwise correlations between gene expression levels were tested using the nonparametric Spearman test. We used the commercial software Genomatix Pathway System [START_REF]Genomatix -NGS Data Analysis & Personalized Medicine[END_REF] to identify significantly enriched pathways and functional themes, among the list of the best correlated genes with genes of interest.

Results

Functional impairment of MD-DCs from SpA patients

To verify that MD-DCs were differentiated from comparable monocyte populations between SpA patients and controls, we analyzed monocytes subsets in both groups. As shown in an additional file, the distribution of the classical, intermediate and non-classical circulating subsets among CD45 + PBMCs was similar in patients and controls (see Additional file 4: Figure S1).

The cells obtained after 7 days of differentiation in culture in patients and controls were uniformly CD14 -, CD11c + , CD40 + , HLA-DR + , CD83 dim and CD86 + , consistent with a MD-DC phenotype. Moreover, HLA-DR, CD80 and CD86 expression increased and CD83 was induced on MD-DCs after LPS stimulation for 24 hours (Additional file 5: Figure S2).

In the rat model of SpA, HLA-B27 transgenic DCs have a decreased capacity to stimulate allogeneic T cells. Thus, we tested the capacity of MD-DCs from HLA-B27 + SpA patients to prime allogeneic CD4 + T cells using an in vitro proliferation assay. A weaker proliferation of CD4 + T cells was observed with SpA MD-DCs as compared to controls, which was statistically significant before and after 6 hours of LPS exposure (P <0.05 and P <0.01, respectively; Figure 1).

Identification of differentially expressed genes in MD-DCs between patients and controls

To investigate the mechanisms underlying the functional defect of MD-DCs from SpA patients, we then compared gene expression levels in those cells between SpA and healthy donors upon stimulation with LPS. After filtering on gene expression level and gene annotation, 13,021 genes (57.3% of the genes present in the microarray) were kept for further analysis. Expression data for MD-DCs either left unstimulated or stimulated with LPS for 6 or 24 hours were first globally compared between SpA patients and healthy subjects. Unsupervised hierarchical clustering showed that samples were perfectly grouped according to LPS stimulation time points (Additional file 6: Figure S3) but not to the disease status. Linear modeling nevertheless identified 81 genes differentially expressed between patients and controls at any time point, with a nominal P-value <0.01 and a fold change below 0.66 or greater than 1.5 (Additional file 7: Table S4). Of this set of genes, 61 were downregulated and 20 were upregulated in patients.

Using the extraction of the dataset corresponding to these 81 genes, LPS stimulation time points were perfectly separated by unsupervised hierarchical clustering. Furthermore, patients and controls were now clearly discriminated, with only four misclassifications (one at H0, two at H6 and one at H24) (Figure 2).

Four candidate genes were chosen for validation study on the basis of their P-value, fold change and biological relevance: ADAMTS15, CITED2, F13A1 and SELL. For all four genes, the qRT-PCR data produced with the RNA samples used for the microarray study confirmed the significant differences in expression and the direction of changes between SpA patients and controls (Table 1, Figure 3): ADAMTS15, F13A1 and SELL were significantly upregulated in SpA samples, whereas CITED2 was downregulated.

Co-regulated genes in SpA patients

Examination of gene expression data indicated that some of them displayed similarity in their expression profile, suggesting co-regulation. Indeed, we found a highly significant negative correlation between ADAMTS15 and CITED2 expression levels at H6 in both patient and control groups (combined data: rho = -0.75; P = 0.0003; Figure 4).

Considering the transcriptional co-activator function of CITED2 and its downregulation in SpA MD-DCs, we sought whether genes other than ADAMTS15 shared a similarly co-regulated expression pattern in patient MD-DCs. We tested the correlation of the genome-wide expression matrix of patient MD-DCs with CITED2 expression in a pairwise fashion. A set of 222 genes was found significantly correlated with CITED2 (P <0.01). They were analyzed with the Genomatix web platform in order to interrogate their biological relevance. This revealed a significant enrichment of genes belonging to the Wnt signaling pathway (P = 2.48 × 10 -4 ) in the SpA group. In the controls, there was no significant correlation between expression levels of CITED2 and the four genes identified in the Wnt pathway (that is, WNT1, WNT10B, FZD4 and ROR2) (Figure 5). This co-regulation pattern involving four genes of the Wnt signaling pathway was therefore specifically associated with the disease.

Discussion

Starting from the hypothesis that aberrant DCs function could play a critical role in the development of SpA, we showed here for the first time that the capacity of MD-DCs to stimulate allogeneic CD4 + T cell response was impaired in SpA patients compared to controls. The genome-wide transcriptome of these MD-DCs elicited by LPS stimulation revealed 81 genes differentially expressed, 4 of which have been validated by RT-qPCR. Finally, we identified the Wnt signaling pathway as dysregulated in patients.

Our choice to study in vitro-differentiated DCs rather than DCs purified from peripheral blood had the advantage to remove those cells from in vivo influences, such as non-specific acute inflammation and drug therapy. Therefore, the differences that we observed between DCs from SpA and controls were more likely to be intrinsic to the DC and could be involved in the disease mechanism.

First, we demonstrated that MD-DCs differentiated from similar monocyte populations in both groups. Then, we observed that the ability of SpA MD-DCs to stimulate allogeneic CD4 + T cell proliferation was markedly defective, as compared to controls, as previously shown in SpAprone HLA-B27 transgenic rats. The molecular basis for this impaired DC function has yet to be elucidated. To investigate this impairment, we have studied the time-dependent global transcriptome of MD-DCs in response to LPS stimulation. A number of microarray-based studies have previously been undertaken in AS, on whole blood cells [START_REF] Gu | Identification of RGS1 as a candidate biomarker for undifferentiated spondylarthritis by genome-wide expression profiling and real-time polymerase chain reaction[END_REF][START_REF] Duan | Gene expression profiling reveals a downregulation in immune-associated genes in patients with AS[END_REF][START_REF] Assassi | Whole-blood gene expression profiling in ankylosing spondylitis shows upregulation of toll-like receptor 4 and 5[END_REF], PBMCs [START_REF] Pimentel-Santos | Whole blood transcriptional profiling in ankylosing spondylitis identifies novel candidate genes that might contribute to the inflammatory and tissue-destructive disease aspects[END_REF][START_REF] Sharma | Insights in to the pathogenesis of axial spondyloarthropathy based on gene expression profiles[END_REF] or macrophages [START_REF] Smith | Gene expression analysis of macrophages derived from ankylosing spondylitis patients reveals interferon-gamma dysregulation[END_REF]. However, no transcriptomic analysis had been conducted on DCs. Of note, as we worked on a homogeneous cell population rather than on cell mixture, our results were not confounded by variations in the composition of the cell populations [START_REF] Plessy | Population transcriptomics with single-cell resolution: a new field made possible by microfluidics: a technology for high throughput transcript counting and data-driven definition of cell types[END_REF]. The time-dependency of our gene profiling is also a novel aspect, rarely investigated in humans.

Our analysis revealed 81 genes differentially expressed in resting and/or LPS-stimulated MD-DCs between SpA patients and controls. Of note, we did not identify a reverse interferon signature, such as in HLA-B27 transgenic rat splenic DCs [START_REF] Fert | Reverse interferon signature is characteristic of antigen-presenting cells in human and rat spondyloarthritis[END_REF] or in SpA patients monocytederived macrophages [START_REF] Smith | Gene expression analysis of macrophages derived from ankylosing spondylitis patients reveals interferon-gamma dysregulation[END_REF]. Such a discrepancy could be explained by differences in culture conditions. Here, we used IL-4 to differentiate MD-DCs, a cytokine that opposes interferon and interferon-induced gene expression and this might have blunted interferon-related differences between groups. Based on the magnitude of their variation between patients and controls, the level of statistical significance, and their biological relevance, we selected four of these genes that we validated using qRT-PCR.

We found an increased expression of ADAMTS15 (A Disintegrin And Metalloproteinase with ThromboSpondin motif) in SpA. Although the ADAMTS family has not yet been implicated in this disorder, numerous studies have identified a role for metalloproteinases (MMPs) in SpA susceptibility and severity [START_REF] Zhu | Matrix metalloproteinase expression in the spondyloarthropathies[END_REF][START_REF] Maksymowych | Serum matrix metalloproteinase 3 is an independent predictor of structural damage progression in patients with ankylosing spondylitis[END_REF][START_REF] Mattey | Association of cytokine and matrix metalloproteinase profiles with disease activity and function in ankylosing spondylitis[END_REF]. Of note, ADAMTS15 was shown to be expressed in the joint with decreased expression in osteoarthritis [START_REF] Kevorkian | Expression profiling of metalloproteinases and their inhibitors in cartilage[END_REF]. Thus, ADAMTS15 could be implicated in the cartilage and/or bone turnover that takes place during joint inflammation, such as in SpA.

CITED2 functions as a context-dependent transcriptional modulator to up-or downregulate the expression of specific genes [START_REF] Bragança | Human CREB-binding protein/p300-interacting transactivator with ED-rich tail (CITED) 4, a new member of the CITED family, functions as a co-activator for transcription factor AP-2[END_REF][START_REF] Lou | Negative feedback regulation of NF-κB action by CITED2 in the nucleus[END_REF][START_REF] Goodman | CBP/p300 in cell growth, transformation, and development[END_REF]. Here, the expression of CITED2 was downregulated in SpA MD-DCs. We further showed an inverse correlation between CITED2 and ADAMTS15 expression after 6 hours of LPS treatment. Interestingly, similar inverse correlation has previously been reported between CITED2 and several MMP family members [START_REF] Bai | A role for CITED2, a CBP/p300 interacting protein, in colon cancer cell invasion[END_REF][START_REF] Lee | Identification of CITED2 as a negative regulator of fracture healing[END_REF]. MMPs constitute a very important group of proteolytic enzymes in joint tissues. Thus, it has been suggested that CITED2 exerted chondroprotective effects through MMP downregulation [START_REF] Leong | Physiological loading of joints prevents cartilage degradation through CITED2[END_REF]. The downregulation of CITED2 and the inverse correlation between CITED2 and ADAMTS15 expression that we observed here suggest that DCs or other phagocytic mononuclear cells, such as macrophages and osteoclasts, could be implicated in SpA joint resorption through heightened metallopeptidase activity.

In silico pathway analysis conducted on CITED2 coexpressed genes highlighted downregulation of several factors belonging to the canonical (that is, WNT1, WNT10B and FRZL4) and non-canonical (ROR2) Wnt signaling pathways, of potential relevance for SpA pathogenesis. First, factors of the canonical Wnt pathway were shown to play a crucial role either in bone formation or destruction in inflammatory arthritis [START_REF] Goldring | Eating bone or adding it: the Wnt pathway decides[END_REF][START_REF] Laine | WNT1 mutations in early-onset osteoporosis and osteogenesis imperfecta[END_REF], such as WNT10b [START_REF] Bennett | Regulation of osteoblastogenesis and bone mass by Wnt10b[END_REF][START_REF] Franck | Evaluation of bone mineral density, hormones, biochemical markers of bone metabolism, and osteoprotegerin serum levels in patients with ankylosing spondylitis[END_REF]. On the other hand, activation of canonical Wntβ-catenin signaling in DCs was shown to concur to regulatory T cell differentiation and conversely to inhibit Th17 differentiation [START_REF] Manicassamy | Activation of beta-catenin in dendritic cells regulates immunity versus tolerance in the intestine[END_REF]. Moreover, ROR2 is a signaling component of the non-canonical Wnt pathway acting on actin cytoskeleton to stimulate cell migration [START_REF] Nishita | Filopodia formation mediated by receptor tyrosine kinase Ror2 is required for Wnt5a-induced cell migration[END_REF][START_REF] Nomachi | Receptor tyrosine kinase Ror2 mediates Wnt5a-induced polarized cell migration by activating c-Jun N-terminal kinase via actin-binding protein filamin A[END_REF]. Interestingly, altered T cell stimulation was linked to defective cytoskeleton dynamics in HLA-B27 transgenic rat DCs [START_REF] Dhaenens | Dendritic cells from spondylarthritis-prone HLA-B27-transgenic rats display altered cytoskeletal dynamics, class II major histocompatibility complex expression, and viability[END_REF].

F13A1 encodes the coagulation factor XIII A subunit, a transglutaminase enzyme. It has multiple extra-and intracellular functions, including a role in cartilage and bone development. Hence, increased factor XIII A expression has previously been associated with cartilage ageing and degenerescence [START_REF] Johnson | Interleukin-1 induces pro-mineralizing activity of cartilage tissue transglutaminase and factor XIIIa[END_REF]. Fibrin crosslinking by factor XIII is of crucial importance not only for hemostasis, but also for inflammation. For instance, factor XIII A-subunit genotype was shown to influence C-reactive protein levels during inflammation in rheumatoid arthritis (RA) [START_REF] Hoppe | Fibrinogen and factor XIII A-subunit genotypes interactively influence C-reactive protein levels during inflammation[END_REF].

Finally, SELL encodes the lymphocyte homing receptor L-selectin/CD62L, one of the major adhesion molecules, which regulates entry of neutrophils and monocytes into inflamed tissues and contributes to the severity of joint inflammation in experimental arthritis [START_REF] Sarraj | Expression of CD44 and L-selectin in the innate immune system is required for severe joint inflammation in the proteoglycan-induced murine model of rheumatoid arthritis[END_REF]. Thus, upregulation of SELL, as shown here in DCs could well participate in joint inflammation in SpA. Interestingly, this gene was also found to be upregulated in splenic DCs from the HLA-B27/human β2-microglobulin transgenic rat [START_REF] Fert | Reverse interferon signature is characteristic of antigen-presenting cells in human and rat spondyloarthritis[END_REF].

Conclusions

Results of the present study reveal the defective functional capacity of DCs from SpA patients, as compared to controls. Furthermore, our results demonstrate significant changes in MD-DCs gene expression upon LPS stimulation that may be inherent to SpA patients. Some of them (that is, decreased Wnt signaling) could account for the altered DC function that provided a rationale for the present study. Others (that is, upregulation of ADAMTS15, F13A1 and SELL) would concur to reinforce tissue inflammation and/or damage. Finally, the co-expression of CITED2 co-transcriptional factor with several of the foregoing genes (that is, the Wnt signaling pathway and ADAMTS15) supports the hypothesis that a coordinated deregulation taking place in DCs may play an important role in SpA pathogenesis.

Figure 1

 1 Figure 1 Impaired allogeneic T cell stimulatory capacity of monocyte-derived dendritic cells (MD-DCs) from spondyloarthritis (SpA) patients. MD-DCs from SpA patients and healthy controls that were left unstimulated at baseline (H0), or were stimulated with lipopolysaccharide (LPS) for 6 h or 24 h (H6 and H24, respectively) were tested for their capacity to stimulate in vitro allogeneic CD4 + T cells. Results are expressed as 3 H-deoxythymidine incorporation after 6 days of mixed lymphocyte reaction, in counts per minute (CPM) ratio. Bars represent the mean proliferation index and standard error of the mean induced by MD-DCs from 19 SpA patients and 24 healthy controls. The Wilcoxon test showed significant differences at H0 (*P <0.05) and H6 (**P <0.01). MLR, mixed lymphocyte reaction.

Figure 2

 2 Figure 2 Unsupervised hierarchical clustering of 57 samples based on microarray expression levels of the 81 differentially expressed genes between spondyloarthritis (SpA) patients and controls. Each time point is color-coded (baseline (H0): blue, 6 h (H6): red, 24 h (H24): green). Samples are clustered on the horizontal axis (top row: SpA patients, bottom row: healthy controls), with the length on the vertical axis representing the degree of correlation between samples.

Figure 3

 3 Figure 3 Scatter plots showing qRT-PCR expression levels in spondyloarthritis (SpA) and controls of the 4 selected genes: ADAMTS15 (A), CITED2 (B), F13A1 (C) and SELL (D). The x-axis of the plots represents the three lipopolysaccharide-stimulation time points (baseline (H0), 6 h (H6) and 24 h (H24)) and the y-axis shows the log2 of gene expression level normalized with housekeeping genes (ΔΔCt). *P <0.05, ***P <0.0005, for SpA versus controls.

Figure 4

 4 Figure 4 Correlation of qRT-PCR expression levels between CITED2 and ADAMTS15. Expression is shown for H6 time-point in mixed spondyloarthritis (open circles) and control (black triangles) samples. P-value (P) and nonparametric Spearman correlation coefficient value (rho) are shown.

Figure 5

 5 Figure 5 Co-expression of the Wnt signaling pathway with CITED2 in spondyloarthritis (SpA). Graphical view of the Wnt pathway identified by studying genes co-expressed at 6 h (H6) with CITED2 in the whole microarray dataset in patients (A) and in controls (B). Numbers and color code (scale on the right side) indicate nonparametric Spearman correlation coefficient of the network's gene with CITED2.

  2010 GEMISA) and the financial support of the Investissements d' Avenir programme ANR-11-IDEX-0005-02, Sorbonne Paris Cite, Laboratoire d'excellence INFLAMEX. Alice Talpin was supported by a grant from Arthritis Fondation Courtin. Félicie Costantino was supported by a grant from the Société Française de Rhumatologie (SFR). Charles-Antoine Dutertre was supported by the Agence Nationale de Recherche sur le SIDA et les Hépatites Virales (ANRS). Sonia Amraoui had a fellowship from Paris-Diderot University.

Table 1 qRT

 1 -PCR expression analysis of four selected genes in monocyte-derived dendritic cells from spondyloarthritis patients, as compared to controls -value is reported considering the global effect, the effect of the disease status only or the effect of the disease status adjusted for time. Bold text indicates significant P-values and at what time the modulation of gene expression was significantly different between patients and controls.

	Gene	Microarray qRT-PCR
		P-value	P-value		Fold change
			Disease Disease/time H0	H6	H24
	ADAMTS15 0.0003	0.0005 0.04	2.2	5.63 2.2
	CITED2	0.0007	0.21	0.031	0.94 0.51 0.99
	F13A1	0.00006	0.019	0.001	1.32 1.25 5.05
	SELL	0.006	0.012	0.8	3.57 3.47 2.97

P

http://arthritis-research.com/content/16/4/417

Acknowledgements

We gratefully acknowledge the contribution of Pr A Hosmalin (Director of Antigen presentation by Dendritic cell team), for fruitful scientific discussions and providing antibodies for monocyte subset characterization. This work was supported by a grant from Agence Nationale de la Recherche (grant ANR

Additional files

Additional file 1: Table S1. Characteristics of the study patients.*SpA, spondyloarthritis; BASDAI, Bath Ankylosing Spondylitis Disease Activity Index; NSAID, non-steroidal anti-inflammatory drug; TNF, tumor necrosis factor. ND, not done. *The registered manifestations correspond to those present at the time of examination, or retrieved from past-medical history. **Refers to radiographic sacroiliitis ≥ grade II bilateral or grade III unilateral. ***Six patients are common to both studies. ****Data available for nine patients.

Additional file 2: Table S2. Characteristics of the study healthy controls.

Additional file 3: Table S3. Nucleotide sequence of the PCR primers.

Additional file 4: Figure S1. Comparison of monocyte subsets distribution among peripheral blood mononuclear cells (PBMCs) between spondyloarthritis (SpA) patients and healthy controls (Ctrl) by six-color flow cytometry. Distribution of monocyte subsets among CD45 + PBMC was studied in five healthy donors (clear boxes) and 10 HLA-B27 + SpA patients (gray boxes). Results are represented as boxes, bars indicate medians. The mean age of healthy donors was 42 years at the time of the study and 40% of them were men. Characteristics of the patients are shown in Additional file 1: Table S1.

Additional file 5: Figure S2. Phenotypic characterization of purified monocytes and monocyte-derived dendritic cells (MD-DCs) by flow cytometry. Results of one healthy donor representative of study subjects are shown. Monocytes (D0) were CD14 + , CD11c + , CD40 -, HLA-I + , HLA-DR + , CD80 -, CD83 -and CD86 -. MD-DCs (D7) were CD14 -, CD11c + , CD40 + , HLA-I + , HLA-DR + , CD80 -, CD83 dim , and CD86 + . HLA-DR, CD83 and CD86 expression increased, and CD80 was induced on MD-DCs after lipolpolysaccharide (LPS) stimulation for 24 h (H24). Clear plot represents control isotype and gray plot the tested antibody.

Additional file 6: Figure S3. Unsupervised hierarchical clustering of 57 samples based on whole-genome gene expression levels. Each time point is represented by a color (baseline (H0): blue, 6 h (H6): red, 24 h (H24): green). Samples are clustered on the horizontal axis (top row: spondyloarthritis (SpA) patients, bottom row: healthy controls).

Additional file 7: Table S4. List of the genes differentially expressed in monocyte-derived dendritic cells (MD-DCs) between spondyloarthritis (SpA) and controls, ranked by fold change.**The criteria for the inclusion of genes in this table are described in Patients and Methods. Global values are LInear models for microarray data (LIMMA) values. The four genes selected for qRT-PCR validation are in bold text.

Abbreviations AS: ankylosing spondylitis; CPM: counts per minute; DC: dendritic cell; EBV: Epstein-Barr virus; GM-CSF: granulocyte-macrophage colony-stimulating factor; IL: interleukin; LIMMA: LInear models for microarray data; LPS: lipopolysaccharide; mAb: monoclonal antibody; MD-DC: monocyte-derived dendritic cell; MLR: mixed lymphocyte reaction; MMP: metalloproteinase; PBMC: peripheral blood mononuclear cell; qRT-PCR: real-time quantitative PCR; SpA: spondyloarthritis; TNF: tumor necrosis factor.

Competing interests

The authors declare that they have no competing interests.

Authors' contributions

All authors were involved in drafting the article or revising it critically for important intellectual content, and all authors read and approved the final manuscript. GC had full access to all of the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. Study conception and design: AT, FC, HJG, MB, and GC. Acquisition of data: AT, FC, NB, AL, SJ, and SA. Analysis and interpretation of data: AT, FC, NB, FL, FD, CAD, HJG, MB, and GC. The corresponding author confirms that all the individuals listed as authors fulfil the uniform authorship credit requirements for manuscripts submitted to medical journals, that is, that they all contributed to the manuscript based on (1) substantial contributions to conception and design, acquisition of data, or analysis and interpretation of data; (2) drafting the article or revising it critically for important intellectual content; and (3) final approval of the version to be published.