nD variational restoration of curvilinear structures with prior-based directional regularization Odyssée Merveille, Benoît Naegel, Hugues Talbot, Nicolas Passat

Abstract-Curvilinear structure restoration in image processing procedures is a difficult task, which can be compounded when these structures are thin, i.e. when their smallest dimension is close to the resolution of the sensor. Many recent restoration methods involve considering a local gradient-based regularization term as prior, assuming gradient sparsity. An isotropic gradient operator is typically not suitable for thin curvilinear structures, since gradients are not sparse for these. In this article, we propose a mixed gradient operator that combines a standard gradient in the isotropic image regions, and a directional gradient in the regions where specific orientations are likely. In particular, such information can be provided by curvilinear structure detectors (e.g. RORPO or Frangi filters). Our proposed mixed gradient operator, that can be viewed as a companion tool of such detectors, is proposed in a discrete framework and its formulation / computation holds in any dimension; in other words, it is valid in Z n , n ≥ 1. We show how this mixed gradient can be used to construct image priors that take edge orientation as well as intensity into account, and then involved in various image processing tasks while preserving curvilinear structures. Experiments carried out on 2D, 3D, real and synthetic images illustrate the relevance of the proposed gradient and its use in variational frameworks for both denoising and segmentation tasks.

Index Terms-directional gradient, proximity operator, nD images, variational framework, filtering, segmentation

I. INTRODUCTION

C URVILINEAR structures are among the most difficult to preserve when carrying out image processing tasks. By curvilinear, we mean objects with a spatial dimension d that can be considered strictly lower than the dimension n of the space in which they are embedded. In this article, we are interested in objects of dimension d = 1 within spaces of dimension n ≥ 2. For n = 2, such objects can be, for instance, roads in remote sensing images, vessels in eye fundus imaging; for n = 3, they can be fibres in composite materials, complex vascular networks in medical images; for n = 4, they can be tracked trajectories of moving objects in 3D+time images; etc.

The preservation of curvilinear structures is challenging for various reasons. They are generally sparsely distributed within images, due to their low dimension. They are also often thin structures, with a thickness similar to the image resolution. In addition, they can have complex geometry and topology, with high curvatures, tortuosity, junctions and bifurcations, etc. All these properties make observing curvilinear structures difficult.

Odyssée Merveille and Benoît Naegel are with Université de Strasbourg, CNRS, ICube, Strasbourg, France ({merveille,b.naegel}@unistra.fr).

Hugues Talbot is with Centralesupélec, CVN, Paris, France (hugues.talbot@centralesupelec.fr).

Nicolas Passat is with Université de Reims Champagne-Ardenne, CReSTIC, Reims, France (nicolas.passat@univ-reims.fr).

In particular, they are easily corrupted by noise, and they often suffer from partial volume effect, due to the above mentioned properties. In this context, many existing image processing methods cannot efficiently discriminate them from noise and artifacts.

For tackling these issues, a relevant solution consists of guiding the process by information related to the location and geometric properties (size, orientation) of curvilinear structures. The robustness of the framework then depends on [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF] the reliability of this prior information, and (2) the ability to efficiently embed and use this information for improving the restoration or segmentation process.

Regarding point (1), several contributions have been specifically devoted to compute information of curvilinear structures from nD images (often with n = 2 or 3). Two representative examples of such contributions are RORPO [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF] and Frangi vesselness [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], that both allow to get access to orientation information on putative curvilinear structures. A complete state of the art is beyond the scope of this article; a brief survey is proposed in Section II-A.

In this article, our purpose is related to point [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF]. Indeed, we aim at developing a framework that explicitly models and uses information about curvilinear structures for improving their efficient processing. This framework is indeed a companion tool of the curvilinear structure estimators of point [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF]. Our approach relies on a variational paradigm, which is versatile in terms of applications and image dimensions.

Many variational formulations stem from a Maximum a Posteriori Bayesian interpretation that is expressed in a sum of two terms, namely a data-fidelity term and a regularization term. The regularization term corresponds to an image model. Among these, Total-Variation-like approaches aim at regularizing the result versus image noise, by minimizing the overall gradient of the segmentation result. Resulting from several decades of research, this classical regularization term usually provides good results except on curvilinear structures. The reasons of this failure are discussed in Section III-A.

In this context, we define a mixed gradient, which merges a standard gradient and a directional gradient that derives from the directional information provided by any curvilinear structure estimator, for instance [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF]. This mixed gradient is then embedded in a variational formulation to form a directional regularization term adapted to curvilinear structures. In addition, we propose an algebraic formulation of this directional regularization, providing a unified definition, independent from the image dimension.

This article is an extended and improved version of the conference paper [START_REF] Merveille | A variational model for thin structure segmentation based on a directional regularization[END_REF]. It is organized as follows. In Section II, we propose a brief survey of curvilinear structure processing. In Section III, we describe the drawbacks of standard gradientbased regularization, and we summarize our strategy for coping with the identified issues. In Section IV, we propose our gradient definitions. In Section V, we show how both standard and directional gradients can be expressed in a matrix formalism that remains homogeneous in any dimensions. In Section VI, experiments involving our gradient-based regularization are proposed for various kinds of 2D / 3D and synthetic / real images. Concluding remarks are provided in Section VII.

II. CURVILINEAR STRUCTURE PROCESSING

The state of the art of curvilinear structure processing and analysis is, amusingly, both thin and wide. It is thin because it is not a mainstream topic in image processing, and it is wide because applications cover many different fields. In particular vascular imaging -i.e., 2D eye fundus imaging and 3D, 3D+time angiographic imaging-is certainly the object of the most intensive activity in this context. The bibliography published over the 25 last years took advantage of most paradigms classically proposed in image or signal processing and, more recently, in machine learning.

The short bibliographic discussion proposed hereafter is necessarily incomplete and focuses on the topics that are the most related to our proposed contributions. The readers interested in a more general discussion will find a global survey in [4, Chapters 1, 2]; a medical-oriented survey may also be found in [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF], [START_REF] Moccia | Blood vessel segmentation algorithms -Review of methods, datasets and evaluation metrics[END_REF].

A. Curvilinear structure detection / description

Before actually processing / analysing a curvilinear structure, it is often relevant to make it more easily detectable in the image. This can be viewed as a filtering task, by improving the signal-to-noise and / or signal-to-background ratio. It can also be viewed as a more semantic task, by determining higher level information, for instance, local size and / or orientation of the curvilinear structure. In this context, two main families of approaches have been developed.

The first relies on linear operators, based on local, differential analysis of images. In particular, the analysis of secondorder derivatives of 3D images were proposed in [START_REF] Sato | 3D multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF], [START_REF] Lorenz | Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2D and 3D medical images[END_REF]. In these pioneering works, the eigenvectors of multiscale Hessian matrices and their associated eigenvalues are analysed to characterise blobs (3D), planar (2D) and curvilinear (1D) structures as well as their scale and orientation. This led to the proposal of measures combining differential information into heuristic formulations. The measure proposed in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] is often considered the current gold-standard. Several variants have been proposed since then, for instance in [START_REF] Krissian | Model-based detection of tubular structures in 3D images[END_REF], [START_REF] Danielsson | Efficient detection of seconddegree variations in 2D and 3D images[END_REF]. Alternatively, steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF] can be expressed, for similar purposes, in terms of a linear combination of basis filters. As such, they are often used to detect oriented features such as curvilinear structures. A framework for 3D steerable filters was, in particular, proposed in [START_REF] Derpanis | Three-dimensional nth derivative of Gaussian separable steerable filters[END_REF], using a n th Gaussian derivative basis filter.

Anisotropic diffusion methods [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF] were also proposed to filter curvilinear structures. Of particular interest is the Coherence Enhancing Diffusion (CED) filter [START_REF]Coherence-enhancing diffusion filtering[END_REF], designed to perform actual anisotropic diffusion by adapting the diffusion direction using the structure tensor, i.e the tensor product of the gradient. More recently, the Hybrid Diffusion using Continuous Switch (HDCS) [START_REF] Mendrik | Noise reduction in computed tomography scans using 3-D anisotropic hybrid diffusion with continuous switch[END_REF] proposed to combine CED and a diffusion filter with actual anisotropic effects, namely Edge-Enhancing Diffusion (EED) [START_REF] Weickert | Anisotropic diffusion in image processing[END_REF]. This combination allows for a strong denoising effect both in the background and inside curvilinear structures, while still preserving their edges. The direction estimation used by both methods depends on the structure tensor which is a reliable estimator for large scale curvilinear structures, but is less robust at small scales. In [START_REF] Welk | Locally analytic schemes: A link between diffusion filtering and wavelet shrinkage[END_REF], authors consider local analytical and semi-analytical schemes that can be used for implementing true anisotropic diffusion filters such as EED. They also make an interesting link with Haar wavelets shrinkage, allowing for a multiscale approach.

The second family relies on nonlinear approaches. In particular, notions of optimal path detection in graph and mathematical morphology were involved in the development of these approaches. At the frontier between these two domains, geodesic paths [START_REF] Cohen | Segmentation of 3D tubular objects with adaptive front propagation and minimal tree extraction for 3D medical imaging[END_REF] were introduced to consider long-range, non-local interactions while still coping with the constraints of thin objects, in particular noise. A curvilinear object detector was also proposed in [START_REF] Rouchdy | Geodesic voting for the automatic extraction of tree structures. Methods and applications[END_REF] using geodesic voting, similar to path density. In [START_REF] Vincent | Minimal path algorithms for the robust detection of linear features in gray images[END_REF], a notion of local optimal path was pioneered. Its purpose is to restrict the research to a given distance, and in a given cone of orientations. This paradigm led to the development of a notion of path opening [START_REF] Heijmans | Path openings and closings[END_REF], enabling a higher flexibility in geometry and size, while preserving a 1D semantics. Algorithmic efforts were conducted to make such approach computationally efficient [START_REF] Hendriks | Constrained and dimensionality-independent path openings[END_REF] and robust to noise [START_REF] Talbot | Efficient complete and incomplete path openings and closings[END_REF], [START_REF] Cokelaer | Efficient robust d-dimensional path operators[END_REF], leading to a notion of robust path opening. In [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF], [START_REF] Merveille | 2D filtering of curvilinear structures by ranking the orientation responses of path operators (RORPO)[END_REF], the notion of RORPO (ranking the orientation responses of path operators), finally built upon these notions, in order to provide a semi-global, nonlinear alternative to the Hessian based approaches in the 2D and 3D cases.

B. Variational paradigm

In this article, we consider the classical variational image restoration problem expressed as the minimization of a twoterm energy defined as follows:

f = arg min f E data (f ) + λE reg (f ) (1)
where f is the restored image, E data is the data fidelity term and E reg is the regularization term. This framework was popularized in the 1990s with active contour, level sets and image restoration models. Mathematically speaking, image restoration is an ill-posed inverse problem. For solving it, it is necessary to impose some regularity on the solution. In [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF], a quadratic regularization term is used. This is highly efficient, however, this can generate blurring effects. This quadratic regularization can be replaced by a 1 gradient norm, called total variation (TV), which better preserves edges [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF].

This framework is very flexible and can be used with sparse and/or a blurring operator in the data fidelity for non-blind deconvolution, for instance in the context of alpha-matting [START_REF] Rhemann | A spatially varying PSF-based prior for alpha matting[END_REF] and super-resolution [START_REF] Dong | Image deblurring and superresolution by adaptive sparse domain selection and adaptive regularization[END_REF], [START_REF] Yang | Image super-resolution via sparse representation[END_REF]. With missing image values in the input, it can be used for image inpainting [START_REF] Chan | Simultaneous total variation image inpainting and blind deconvolution[END_REF], [START_REF] Fadili | Inpainting and zooming using sparse representations[END_REF].

By constraining the output image to take a restricted set of values, this framework can be used for segmentation. In particular, the Chan-Vese model [START_REF] Chan | Active active contours without edges[END_REF] divides the image into two regions of piecewise constant intensities. If these two constant values are known, it results in a convex problem that can be solved exactly [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF].

In the context of curvilinear structure segmentation, extensions of the Chan-Vese model were proposed by adding curvilinear priors, for instance, superellipsoids [START_REF] Tyrrell | Robust 3-D modeling of vasculature imagery using superellipsoids[END_REF], B-splines framelet [START_REF] Tai | Wavelet frame based multi-phase image segmentation[END_REF], adaptive dictionaries [START_REF] Rigamonti | Accurate and efficient linear structure segmentation by leveraging ad hoc features with learned filters[END_REF] and elastical regularization [START_REF] El-Zehiry | Vessel segmentation using 3D elastical regularization[END_REF]. In [START_REF] Miraucourt | Variational method combined with Frangi vesselness for tubular object segmentation[END_REF], we also proposed a variational approach for tubular structure restoration. By considering the Frangi measure [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF], we designed an adaptive regularization parameter to avoid intensity loss in curvilinear structures, which is an intrinsic problem of classical variational frameworks. However, this approach can prevent regularization within the curvilinear structures, thus leading to various problems, such as potential disconnections.

In the current manuscript, we propose a framework for defining and embedding a mixed gradient operator coupling the directional and standard gradient. This allows us to more robustly take advantage of prior knowledge related to curvilinear structure analysis. Indeed, regularization can be carried out everywhere in the image while remaining adapted to the geometric context. In addition our proposed variational framework remains fully n-dimensional, thus encompassing any images considered in potential applications.

III. CURVILINEAR STRUCTURE RESTORATION USING DIRECTIONAL REGULARIZATION

In this section, we first explain the limitations of classical regularization terms in the context of curvilinear structure restoration. Then, we give the intuition behind our directional regularization and why it is a better suited term for curvilinear structure restoration. More formal explanations on the directional regularization and its implementation details are exposed in Section IV.

A. Regularization principle and limitations

The regularization term in a variational restoration problem can be interpreted as choosing a solution with desirable properties within a solutions space. A property often desired is to minimize image noise while retaining image content. This property is translated into a regularization term that minimizes the norm of the image gradient, under the assumption that the image content has sparse contours. When this norm is the 1 norm, the regularization term is called total variation (TV):

E reg TV (f ) = ∇f 1 ( 2
)
where ∇f is the discrete gradient of f and • 1 is the 1 norm.

Interpreting this regularization term as a statistical model, it says that gradient intensities follow an exponential distribution. This is indeed observed to be the case for natural images [START_REF] Ruderman | The statistics of natural images[END_REF]. Interpreting it as an image model, it says that natural images should have a sparse gradient, meaning that most variations in images are small, except near a few contours where they can be large. This corresponds to an intuitive piecewise smooth image model. This regularization term is effective at decreasing the image noise in natural images; however it may also suppress image contours when they are not highly contrasted. Moreover, structures with a high perimeter over surface ratio (for n = 2) or a high surface over volume ratio (for n = 3) are highly penalized. For these reasons, classical regularization terms are not a good model for curvilinear structures, which tend to disappear in the resulting restored images.

To cope with this problem, Miraucourt et al. [START_REF] Miraucourt | Variational method combined with Frangi vesselness for tubular object segmentation[END_REF] proposed to include a curvilinear structure position prior in order to regularize more strongly outside the curvilinear structures than inside (see Fig. 1(b)). This strategy effectively prevents curvilinear structures from disappearing, but is not effective at decreasing noise inside the curvilinear structures, where only a weak regularization is applied.

B. Directional regularization motivation

The motivation behind [START_REF] Miraucourt | Variational method combined with Frangi vesselness for tubular object segmentation[END_REF] was to keep the good behavior of the total variation outside curvilinear structures and decrease its effects within them. However, nothing specific was done to regularize and restore the curvilinear structures specifically. We propose to extend the total variation formulation to better preserve curvilinear structures, while keeping its good properties on other structures. To this end, we consider not only a positional prior, as in [START_REF] Miraucourt | Variational method combined with Frangi vesselness for tubular object segmentation[END_REF], but also an orientation prior. Instead of a weaker regularization inside curvilinear structures, we apply a strong regularization but solely along the local curvilinear structure direction. In other words, we propose an intensity and directional spatially variant regularization term (see Fig. 1(c)).

The total variation is an isotropic regularization term as the discrete gradient lacks a privileged direction. Indeed, the discrete gradient is computed by finite differences in all the directions of the Cartesian basis. Our goal is to change the total variation behavior inside curvilinear structures to better preserve them. We define a directional gradient ∇ d f which has a privileged direction, namely the local direction of the curvilinear structure itself.

Based on both isotropic (∇f ) and anisotropic (∇ d f ) gradients, we propose to adapt the total variation formulation to better preserve curvilinear structures, while keeping its good properties on other structures. This new regularization term, called directional regularization is defined by:

E reg (f (x)) = ∇ m f (x) 1 (3)
where:

∇ m f (x) = ∇ d f (x) if x ∈ curvilinear structure ∇f (x) otherwise (4)

C. Implementation

In practice, we aim at solving the following problem (see Eq. ( 1)):

f (x) = arg min f ∈C E data (f (x)) + λE reg (f (x)) = arg min f ∈C E data (f (x)) + λ ∇ m f (x) 1 ( 5
)
where C is a closed convex subset of the image domain, E data (f ) is a convex, differentiable function and ∇ m f 1 is a convex, but non-differentiable function that forbids using a gradient descent approach. For many inverse problems including segmentation, the convex C is a centered unit ball which ensures that output values are consistent with the input image data. It can be shown from [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] that Eq. ( 5) admits a solution given by the proximal point splitting algorithm:

u n+1 = prox γh u n -γ∇E data (u n ) (6)
where γ ∈ (0, +∞) is a step-size parameter and

h = λ ∇ m f (u n ) 1 .
The reader may refer to [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] for a review on proximal splitting point algorithm. Here, we only recall the definition of the proximity operator of a function φ:

prox φ (y) = arg min x φ(x) + 1 2 x -y 2 2 (7)
where • 2 is the 2 norm. There is no generic strategy in the literature, that allows one to compute the proximity operator of any function φ. Nonetheless, several algorithms were proposed to compute the proximity operator of specific functions.

To solve our restoration problem (Eq. ( 5)), we use the iterative method of Eq. ( 6), that requires to compute the proximity operator of γh. In particular, we choose the Beck and Teboulle algorithm [START_REF] Beck | Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems[END_REF], called Fast Gradient Projection (FGP), that was designed to solve this specific proximity operator (see Algorithm 1).

In this algorithm, D is the domain of the image as described below in Eq. ( 8), P 2 is the projection on the 2 unit ball, P C is the projection on the convex set C, γ is set to the Lipschitz constant of the gradient of E data . ∇ and div are the gradient and divergence operators, respectively.

Algorithm 1: FGP algorithm

Data: y ∈ R D , x 0 ∈ R nD , z 1 = x 0 and t 0 = 1 for n ≥ 1 do y n = z n + 1 γλ ∇ P C y -λ div(z n ) x n = P 2 [y n ] t n+1 = 1 + 4t 2 n + 1 2 λ n = t n -1 t n+1 z n+1 = x n + λ n (x n -x n-1 ) return prox γh (y) = P C y -λ div(z N )

IV. GRADIENTS FORMULATIONS

In this section, we define formally the operators presented in the previous section: the gradient, the directional gradient and the mixed gradient.

In the following we define a nD image f (n ∈ N, n ≥ 1) as follows:

f : D → R x → f (x) (8)
where

D = n i=1 [[0, d i -1]] ⊂ Z n and (d i ) i∈[[1,n]]
are the dimensions of the image.

A. Standard gradient

The discrete gradient of an image can be expressed, via finite differences such that:

∇f (x) = f (x + e i ) -f (x) n i=1 (9)
where (e i ) n i=1 is the canonical basis of R n and Z n , namely e i = (δ i,j ) n j=1 for any i ∈ [[1, n]], (with δ i,j the Kronecker delta).

B. Directional gradient

When a curvilinear structure is identified at a point x ∈ D, one may wish to define a gradient operator locally oriented in the direction of this structure, in order to allow for noise removal without altering this structure.

Let d(x) be the unit vector lying in the direction of the curvilinear structure observed at x. We assume that d(x) is oriented in the half-space of R n such that d n ≥ 0 (see discussion in Appendix A). We define the directional gradient ∇ d f : D → R n via finite differences such that:

∇ d f (x) = f (x + d(x)) -f (x) .d(x) (10)
However, f (x + d(x)) is generally undefined, since d(x) ∈ R n while f is a function on D ⊂ Z n . It is then necessary to consider an interpolation of f on the part of R n associated to D. In particular, the standard n-linear interpolation can be considered, namely: where: 2 shows two examples of interpolation in the case n = 2.

f (x + d(x)) = u∈[[-1,1]] n-1 ×[[0,1]] λ d(x),u f (x + u) (11)
λ d(x),u = n i=1 w(d i (x), u i ) (12)
w(d i (x), u i ) = (1 -|u i |)(1 -|d i (x)|) +|u i |(u i .d i (x) + |d i (x)|)/2 (13)
with x = (x i ) n i=1 , d(x) = (d i (x)) n i=1 , u = (u i ) n i=1 . Fig.

C. Mixed gradient

It is possible to associate to the function f : D → R, a directional vector field d : D → R n that yields a unit vector d(x) for each point x ∈ D if a specific orientation is detected at f (x), and the null vector 0 otherwise.

In particular, the definition of d(x) either as a unit vector providing a specific orientation, or a null vector expressing the absence of such orientation, is conditioned by the curvilinear structure estimator used beforehand. For instance, this can be a thresholding of the RORPO filter or of the Frangi vesselness output. In case of a positive answer, both RORPO and Frangi vesselness provide orientation information.

This vector field provides the information required to compute the directional gradient ∇ d f , see Section IV-B. It is also useful for computing and modeling in a matrix fashion a mixed directional-standard gradient of f .

In particular, let us consider the function α : D → {0, 1} defined, for any x ∈ D by α(x) = d(x) 2 , namely the norm of d(x). This function takes its values in the binary set {0, 1}. It is equal to 1 when a unit vector d(x) is defined at x, and 0 when d(x) is the null vector 0, i.e. when no specific orientation is defined at x.

Then, according to the binary values of α, we can build a mixed gradient ∇ m f from the standard and directional gradients ∇f and ∇ d f as:

∇ m f (x) = ∇f (x) if α(x) = 0 ∇ d f (x) if α(x) = 1 (14)
that is:

∇ m f (x) = α(x)∇ d f (x) + (1 -α(x))∇f (x) (15)
V. MATRIX FORMULATION OF THE DIRECTIONAL REGULARIZATION Even though the vector expression of the mixed gradient presented above (see Eq. ( 15)) is simple, its implementation and manipulation in nD (n > 2) becomes somewhat complex. Indeed, Algorithm 1 requires the definition of our mixed gradient operator and its associated divergence given by the following adjoint relation:

-div p, u S = p, ∇u S n (16)
with S = R D . The higher the dimension of the image, the more complex the divergence definition. This is especially true in terms of limit cases due to the discrete nature of an image.

To tackle these issues, we propose to consider the discrete calculus framework, by proposing a matrix formulation M ∇m of the mixed gradient operator. Thus, the adjoint divergence of this matrix operator, M ∇m , will be simply obtained as the transpose matrix (M ∇m ) T .

A. Vector formulation of functions

A function h : X → Y can be modeled as a vector V h ∈ Y |X| of |X| elements, representing the image of each element of X by h, that is:

V h = [h(x)] x∈X =    h(a) . . . h(z)    (17)
In particular, the function f defined in Eq. ( 8) is modeled by the vector:

V f = [f (x)] x∈D =    f ((0, . . . , 0)) . . . f ((d 1 -1, . . . , d n -1))    =    v f 1 . . . v f D    (18)
where D = |D|.

In particular, for modeling f : D → R as a vector V f ∈ R D , we need a transfer function σ between D ⊂ Z n and [[1, D]] ⊂ Z. Such a function can be defined as:

σ : D → [[1, D]] x = (x i ) n i=1 → 1 + n i=1 ( i-1 j=1 d j )x i (19)
The function σ is bijective. In the sequel, we note

x j = σ -1 (j) ∈ D.
Following the notations of Eq. ( 18), we then have:

f (x) = v f σ(x) (20)
and we can then establish a one-to-one correspondence between the terms of f and V f

B. Matrix formulation of the gradient operator

We are now ready to establish the matrix formulation of the standard gradient operator ∇f of f . We have ∇f : D → R n . Then, it can be defined as a vector V ∇f ∈ (R n ) D of D elements or, equivalently, as a vector V ∇f ∈ R n.D of n.D elements. Practically, this vector can be split into n vectors V ∇f i of D elements, each of them giving the coordinates of the gradient in the subspace R of R n induced by a basis vector e i .

In other words, we have:

V ∇f =    V ∇f 1 . . . V ∇f n    (21)
with, for all i ∈ [[1, n]] (see Eq. ( 9)):

V ∇f i =          ∂f ∂xi (x 1 ) . . . ∂f ∂xi (x j ) . . . ∂f ∂xi (x D )          (22) =         f (x 1 + e i ) -f (x 1 ) . . . f (x j + e i ) -f (x j ) . . . f (x D + e i ) -f (x D )         (23)
Consequently, V ∇f i can be expressed as:

V ∇f i = (G i -I D ) • V f (24)
where I D is the D × D identity matrix and

G i = (g i (j, k)) (j,k)∈[[1,D]] 2
is the D × D matrix defined by:

g i (j, k) = δ k,σ(xj +ei) (25)
For the correct handling of boundary effects on the frontier of the domain D, whenever we have x j + e i / ∈ D we set:

g i (j, k) = δ k,σ(xj ) = δ k,j (26)
In other words, in such a case, we consider a Neumann condition of nullity of the gradient coordinate. (Note that this strategy will be also involved in anisotropy effects correction; see discussion in Appendix A.) Finally, we can then define the vectorial expression of ∇f by combining its n coordinate vectors, as:

V ∇f = M ∇ • V f =                 G 1 . . . G i . . . G n         -         I D . . . I D . . . I D                 • V f (27)
where M ∇ is a n.D × D matrix.

C. Matrix formulation of directional gradient

In a similar way, it is also possible to express a matrix formulation of the directional gradient. This formulation is, however, slightly different from that of the standard gradient. Indeed, for any x ∈ D, ∇f (x) is expressed by its n coordinates with respect to the canonical basis (e i ) n i=1 . In contrast, ∇ d f (x) has only one coordinate with respect to a unit vector d(x) on the 1D line locally oriented in the specific direction observed at x. In particular, this vector d(x) depends on x, and the coordinate system is then spatially variant.

As a consequence, we have to compute only one coordinate of the directional gradient for each x ∈ D, but this computation is more complex than for standard gradient, as it theoretically involves 3 n points of D around x (in practice, only 2 n actually contribute to the result); see Eqs. [START_REF] Freeman | The design and use of steerable filters[END_REF][START_REF] Derpanis | Three-dimensional nth derivative of Gaussian separable steerable filters[END_REF][START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF].

The vectorial representation V ∇ d f of ∇ d f is then expressed as:

V ∇ d f =          ∂f ∂d(x1) (x 1 ) . . . ∂f ∂d(xj ) (x j ) . . . ∂f ∂d(x D ) (x D )          (28) =         f (x 1 + d(x 1 )) -f (x 1 ) . . . f (x j + d(x j )) -f (x j ) . . . f (x D + d(x D )) -f (x D )         (29)
Consequently, V ∇ d f can be written as:

V ∇ d f = M ∇ d • V f = E D • (G d -I D ) • V f (30)
where

I D is the D × D identity matrix, G d = (g d (j, k)) (j,k)∈[[1,n]] 2
is the D × D matrix defined by:

g d (j, k) = u∈[[-1,1]] n δ k,σ(xj +u) .λ d(xj ),u (31)
and

E D = [ε j ] D j=1 ∈ {0, 1} D is defined by: ε j = 1 -δ 0, u∈[[-1,1]] n σ(xj +u) (32)
For the correct handling of boundary effects on the frontier of the domain D, we still extend -by convention-σ out of D by setting σ(x) = 0 for any x ∈ Z n \ D. With this formulation, a useful feature is that the directional gradient ∇ d f vanishes exactly on the border of the domain D, i.e. for points x = (x i ) n i=1 with at least one coordinate x i = 0 or d i -1.

D. Matrix formulation of the mixed gradient

A matrix representation M ∇m of the mixed gradient ∇ m can be formulated from the above matrix representations of ∇ d and ∇, see Eqs. ( 27) and [START_REF] Yang | Image super-resolution via sparse representation[END_REF]. Indeed, we have:

M ∇m =           (1 D -V α ) • (G 1 -I D )
. . .

(1 D -V α ) • (G i -I D ) . . . (1 D -V α ) • (G n -I D ) V α • E D • (G d -I D )           (33)
where 1 D is the vector representation of the constant function x → 1 on D, and V α is the vector representation of the function α. Note that M ∇m is a (n + 1).D × D matrix.

VI. EXPERIMENTS

In this section, we experimentally assess the relevance of the proposed restoration framework. To emphasize the versatility of our approach in terms of dimensionality and application, we first present a segmentation application on 2D images; then we show denoising results on 3D images. For both experiments, the prior information on estimated orientations of putative curvilinear structures was obtained by the RORPO operator [START_REF] Merveille | Curvilinear structure analysis by ranking the orientation responses of path operators[END_REF].

A. Segmentation

The DRIVE database [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF] is a well known dataset of 40 eye fundus images associated with their blood vessels groundtruth. In this section, we apply our restoration framework to the blood vessel segmentation problem. In this case, we use the data fidelity term proposed by Chan et al. [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF], E data = c f , u F , where:

• c f ∈ R D is the function x → (c 1 -f (x)) 2 -(c 2 -f (x)) 2
(with c ∈ R some constant values assumed to correspond to the background and foreground of the image);

• • , • F the Frobenius product. To assess the relevance of our regularization term in a segmentation framework, we compare the classic Chan et al. segmentation model [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] that uses the total variation:

û = arg min u∈[0,1] D c f , u F + λ ∇u 2,1 (34)
with our segmentation model that uses the proposed directional regularization:

û = arg min u∈[0,1] D c f , u F + λ ∇ m u 2,1 (35)
Both standard and directional versions of the method share the same parameters: c 1 , c 2 and λ. We optimized these parameters based on the images of the DRIVE training set.

The optimization scheme and parameters values are presented in Appendix B.

For our experiments, we consider the single-channel version of these images, i.e. the grey-level version of the data. We subtracted the median filter to each image in order to homogenize the image background, since the Chan et al. data fidelity assumes homogeneous background and foreground intensities. An example of such image is shown in Fig. 3(a).

We computed the true positive rate / sensivity (TPR) and true negative rate / specificity (TNR), along with the accuracy 1 (Acc) in the field of view area provided by the DRIVE dataset of each segmentation: 1 In Section VI-B, we will consider the Matthews correlation coefficient (MCC) [START_REF] Matthews | Comparison of the predicted and observed secondary structure of t4 phage lysozyme[END_REF] instead of accuracy, as it is better suited in the case of sparse structures. Here, the choice of accuracy was motivated by the fact that most of the methods benchmarked on the DRIVE dataset were compared using the accuracy measure. [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]. TPR: TRUE POSITIVE RATIO (SENSITIVITY). TNR: TRUE NEGATIVE RATIO (SPECIFICITY). ACC: ACCURACY. TPR TNR Acc Chan [START_REF] Chan | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF] 0.6615 0.9833 0.9423 Our approach 0.6820 0.9809 0.9427 Staal et al. where TP (resp. FP) is the number of true (resp. false) positives and TN (res. FN) is the number of true (resp. false) negatives.

TPR = TP TP + FN TNR = TN TN + FP Acc = TP + TN TP + TN + FP + FN (36)
The mean quantitative results are summarized in Table I and illustrated in Fig. 3. For comparison, we also provide some state of the art results.

Qualitatively, the segmentation improvements mostly concern the distal part of the blood vessels. In particular, our directional regularization successfully reconnects these extremities, as shown in Fig. 3(e-g). This leads to a better connected vessel network, which is a highly desired feature in blood vessel related applications.

Quantitatively, the directional regularization improves the accuracy of the classical Chan et al. segmentation. As the reconnections represent only a few pixels within the image, the accuracy of our directional regularization is only slightly higher than the Chan et al. accuracy even if the improvement is real and significant as can be seen on the images.

Moreover, the accuracy of our method is close to state of the art methods, even though our purpose was not to propose a dedicated segmentation model for retinal images, but only a generic directional regularization term for curvilinear structures.

B. Denoising

In this section, we assess the performance of our method in the case of 3D image denoising. It is difficult to have access to reliable ground-truth for 3D images with curvilinear structures. Thus, we generated synthetic images using the VascuSynth software package2 [START_REF] Hamarneh | VascuSynth: Simulating vascular trees for generating volumetric image data with ground-truth segmentation and tree analysis[END_REF]. We generated 10 images containing curvilinear patterns, organized as a connected acyclic structure, constituting our ground-truth. We then added a Gaussian random field background to these images, in order to simulate undesired, non-homogeneous, smooth, blob-like features. Finally, we introduced a non-additive, spatially variant noise by subjecting each image with 7 different intensities of Poisson noise, determined by the gain g. This resulted in a database of 70 images. An example of such image is presented in Figures 5(b) and 6(b).

We compared our results with three methods: the exact ROF model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], which is a classical variational denoising filter, and is thus used as the baseline experiment showing how our directional regularization improves the curvilinear structures restoration; the hybrid diffuson using continuous switch (HDCS) [START_REF] Mendrik | Noise reduction in computed tomography scans using 3-D anisotropic hybrid diffusion with continuous switch[END_REF], which is a state of the art anisotropic diffusion filter designed to preserve and enhance curvilinear and blob-like structures while filtering the noise; and BM4D [START_REF] Maggioni | Nonlocal transform-domain filter for volumetric data denoising and reconstruction[END_REF], which is the 3D version of BM3D, considered as one of the gold standard denoising algorithms. These experiments were performed using the TubeTK implementation of HDCS 3and the original implementation of BM3D in Matlab 4 .

For each of these methods, we optimized the parameters on the first image of the dataset and for each level of noise. Then, we used this fixed set of parameters on the whole dataset. Parameter optimization was performed by sampling the parameter space for each method in reasonable ranges. The sets of parameters used for each method is provided in Appendix B.

For each method and level of noise, we computed the mean ROC curve over the 10 images as the true positive rate / sensivity (TPR) against the false positives rate (FPR). As the amount of negatives is much higher than the amount of positives in these images, we define the FPR as the number of false positives over the number of positives (P), namely FPR = FP / P. For instance, a FPR of 2 means that there is twice as many false positives as ground truth positives. We also computed the Matthews correlation coefficient (MCC) [START_REF] Matthews | Comparison of the predicted and observed secondary structure of t4 phage lysozyme[END_REF] of the best threshold of each result:

MCC = TP • TN -FP • FN (TP + FP)(TP + FN)(TN + FP)(TN + FN) (37)
The choice of MCC is motivated by the fact that it is well suited for sparse images. The closer to 1 the MCC, the more similar the result to the ground-truth. We chose to analyze these denoising results based on two criteria that are usually considered as segmentation metrics because the ground-truths of this denoising are the binary images of curvilinear structures used for the dataset generation. In this case, the best denoising result is the closest to these binary structures.

A summary of the results is shown in Table II, and the mean ROC curves for a representative noise level (gain g = 15) is shown in Figure 4. From a quantitative point of view, we observe that the MCC scores are significantly higher with our directional regularization, compared to the other three methods irrespective of the noise level. The ROC curves show that our directional regularization outperforms all the compared methods when the FPR is between 0 and 2. For a FPR between 2 and 4, our proposal still outperforms ROF and BM4D, but in this situation HDCS provides the best results.

In other words, our method is the best for preserving the 3D curvilinear structures in the most challenging and useful situations, corresponding to the cases were the TPR is high while the FPR must remains low. This quantitative analysis is qualitatively confirmed by visual inspection of Figures 5 and6. We observe that the reconstruction of the thinnest structures, with thickness close to the image resolution, is better with the directional gradient. This is indeed coherent since the directional regularization is carried out only in orientations that do not lead to alteration of the contours of the structures of interest. BM4D appears to have the best restoration result but it actually oversmooths.

VII. CONCLUSION

In this article, we have proposed a regularization term for curvilinear structure restoration based on the design of a mixed gradient operator. This spatially-variant operator can be computed from orientation measures provided by curvilinear structure detectors. This regularization term has a different behaviour depending on the location in the image. Within curvilinear structures, it regularizes only in the local main direction of the considered structure. Within isotropic structures, it behaves like an isotropic gradient. Using the formalism of discrete calculus, this operator has a consistent behaviour everywhere. In particular it handles border effects gracefully both on the border of the image and at the interface between oriented and isotropic areas. The adjoint divergence operator is immediately derived again thanks to the discrete calculus formulation.

The mixed gradient developed in this work was developed and formalized in a non-dimensional way, thus leading to restoration approaches that can proceed irrespective of the space dimension (2D, 3D, etc.).

We illustrated our approach on segmentation using the classical Chan et al. data fidelity term, and image denoising using the 2 norm fidelity term. However, other formulations may be considered, as long as they include a regularization term expressed as a gradient measure. For instance, it would make sense to use the exact Poisson log-likelihood as data fidelity term for the Poisson noise case [START_REF] Chouzenoux | A convex approach for image restoration with exact poisson-gaussian likelihood[END_REF]. Our next work will augment these strategies with additional priors (for instance considering connectivity), as well as non-local oriented gradient approaches.

Our approach may also be compliant with deep-learning strategies recently investigated in the field of curvilinear structure analysis. Indeed, deep-learning is emerging in this domain, mainly focusing on 2D cases at this stage [START_REF] Litjens | A survey on deep learning in medical image analysis[END_REF], with a specific focus, e.g. on retinal imaging [START_REF] Liskowski | Segmenting retinal blood vessels with deep neural networks[END_REF]. However, very recent works also aim at considering more complex cases, in 3D (bio)medical imaging modalities [START_REF] Tetteh | DeepVesselNet: Vessel segmentation, centerline prediction, and bifurcation detection in 3-D angiographic volumes[END_REF], [START_REF] Haft-Javaherian | Deep convolutional neural networks for segmenting 3D in vivo multiphoton images of vasculature in Alzheimer disease mouse models[END_REF], [START_REF] Sanches | Cerebrovascular network segmentation of MRA images with deep learning[END_REF]. In this context, our mixed gradient may be involved as an additional loss term for strengthening connectivity of thin, curvilinear structures in complex nD images.

As most discrete gradient formulations, the one proposed in this article is anisotropic. Indeed, the standard gradient formulation described in Section IV-A is based on finite differences between values at x and at x+e i for the e i vectors of the canonical basis of Z n (see Eq. ( 9)). The choice of x+e i versus xe i is arbitrary, and leads to compute the gradient in a cone that represents 1/2 n of the neighbourhood of x in Z n . This anisotropy, together with the combined use of two kinds of gradients within a same image, may cause some slight boundary effects at the frontier between regions with / without specific directions. For instance, let us consider two neighbour points x and y such that α(x) = 0 and α(y) = 1 (see Eq. ( 14)), i.e. with a standard gradient at x and a directional gradient at y.

At y, the gradient formulation ∇ m f (y) = ∇ d f (y) will take into account the difference of intensities between I(y) and I(x) only if the direction d(y) is mainly the same as x-y. Then, at y, the global behaviour of ∇ m f is as expected.

By contrast, at x, the gradient formulation ∇ m f (x) = ∇f (x) will take into account the difference of intensities between I(y) and I(x) whenever y-x is one of the e i vectors of the canonical basis. In such case, a gradient value will be computed in x, while it may correspond to the (external) border of a curvilinear structure, where we expect to vanish the gradient. This phenomenon happens, by definition, on one side of the objects, oriented opposite to the canonical vector basis. This may cause, for instance, non-symmetric, blurring effects in segmentation results.

A simple way to get rid of such undesired effects consists of vanishing the i th component of the standard gradient at x whenever x + e i is a point of D where a specific orientation is defined.

Practically, we only have to modify the conditions of validity of Eqs (25-26) as follows. We set:

D 0 = {x ∈ D | α(x) = 0} (38)
D 1 = {x ∈ D | α(x) = 1}
In other words, D 0 (resp. D 1 ) corresponds to the part of D where a standard (resp. directional) gradient is valid (see Eq. ( 14)). Then, we slightly modify the definition of the matrices G i , i ∈ [ [1, n]] (see Eq. ( 24)) constituting the gradient matrix of V ∇f as follows:

g i (j, k) = δ k,σ(xj +ei) if x j + e i ∈ D 0 δ k,j if x j + e i / ∈ D 0 (39)
In other words, we now apply the Neumann condition of nullity of the gradient coordinates on the border of D 0 , i.e. on the frontier of the whole area of non-oriented points, and not only on the boundary of the image. The overall idea remains the same, that is avoiding to compute a classical gradient coordinate when the information carried by one of the two points is not known / not relevant for this kind of gradient.

B. Parameter optimization

In the following, the parameters of each method used in the experiments of Section VI are enumerated, and we present the parameter values used in these experiments.

1) Segmentation: Both compared methods share the same parameters: the background and foreground constants c 1 and c 2 , and the regularization weight λ balancing between the data fidelity and the directional regularization. The constants c 1 and c 2 depend on the initial image. The images were preprocessed to homogenize the background, resulting in a mean background intensity of almost c 1 = 0. We then chose to set the c 2 constant to the Otsu threshold value of each image. Finally we performed a heuristic optimization of the λ parameter by looking for the best value in the range [0.001, 0.5] for each training image for both methods; we chose for each method the most frequent λ value. This resulted in λ = 0.10 for the ROF method and λ = 0.15 for our method. One may note that our directional regularization allows for a higher regularization power without losing performance.

The directional regularization also requires a directional information provided here by the RORPO operator. RORPO has 3 parameters: the minimum scale s, the factor f and the number of scales n. They mostly depend on the dimension of the structure of interest, and the results are robust to parameter variations. We set them to: s = 10, f = 2, n = 5.

2) Denoising: The parameter optimization of this experiment was performed on the first image of the VascuSynth dataset. We sampled the parameter space of each method in reasonable ranges and we kept the set of parameters leading to the best MCC on this image.

• Directional regularization -The directional regularization requires a directional information provided in this work

Fig. 1 .

 1 Fig. 1. Regularization principle. The dark (resp. bright) blue areas represent high (resp. low) isotropic (i.e. with no privileged direction) regularization, whereas the red area represents a directional regularization. (a) The classical regularization is performed in all directions with the same intensity everywhere. (b) The regularization proposed by Miraucourt et al.[START_REF] Miraucourt | Variational method combined with Frangi vesselness for tubular object segmentation[END_REF] acts with a low intensity inside curvilinear structures and a high intensity elsewhere, but always in all directions. (c) Our directional regularization acts with the same intensity inside and outside curvilinear structures; however, inside curvilinear structures, we only regularize along the curvilinear structure direction.

Fig. 2 .

 2 Fig. 2. Illustration of two interpolations (in red and green) in the case n = 2. d A (resp. d B ) is the local orientation of a curvilinear structure at point x (x = x + u 1 ). To compute the directional gradient, the value at point f (x + d A (x)) (resp. f (x + d B (x))) is required.

Fig. 3 .

 3 Fig. 3. Comparison of segmentation results with the classical Chan et al. model (c) and our directional TV model (d) on the first DRIVE image (a) and its ground-truth (b). (e-g) Zoom on distal parts of blood vessels of the Chan model results (bottom) and our proposed model results (top).

Fig. 4 .

 4 Fig. 4. Mean ROC curves over the 10 VascuSynth images with with a Poisson noise of gain g = 15, for FPR values between 0 and 10. (Note that due to the sparsity of the images, the maximum TPR is reached for all methods only for a FPR around 100.)

Fig. 5 .

 5 Fig. 5. Denoising results on a 3D synthetic image (2D slices). (a) Synthetic image generated with VascuSynth [48]. (b) Noisy image generated from (a) with Poisson noise (gain g = 15). (c) Denoising image with total variation (ROF model). (d) Denoising image with directional regularization.

Fig. 6 .

 6 Fig. 6. Denoising results on a 3D synthetic image (maximum intensity projections). (a) Synthetic image generated with VascuSynth [48]. (b) Noisy image generated from (a) with Poisson noise (gain g = 15). (c) Denoising image with total variation (ROF model). (d) Denoising image with directional regularization.

TABLE I QUANTITATIVE

 I SEGMENTATION RESULTS ON THE DRIVE DATABASE

TABLE II QUANTITATIVE

 II DENOISING RESULTS ON VASCUSYNTH IMAGES. THE MEAN MCC VALUE OVER THE 10 IMAGES FOR THE SAME LEVEL OF NOISE IS GIVEN FOR EACH METHOD.

		Noise (g)	Ours ROF [27] HDCS [16] BM3D [50]
		2 0.668	0.659	0.593	0.564
		5 0.580	0.560	0.528	0.460
		10 0.457	0.441	0.438	0.375
		15 0.510	0.442	0.464	0.400
		20 0.463	0.407	0.432	0.379
		30 0.336	0.275	0.321	0.264
		50 0.286	0.227	0.265	0.212
		1.0
		0.8
	TPR = TP/P	0.4 0.6
		0.0 0.2					Directional ROF HDCS BM4D
		0	2	4	FPR = FP/P	6	8	10

TABLE III ROF

 III AND DIRECTIONAL REGULARIZATION PARAMETERS.by the RORPO operator. RORPO has 3 parameters: the minimum scale s, the factor f and the number of scales n. These parameters mostly depend on the size of the structures of interest, but the results are robust to parameter variations. We set them once for all noise levels: s = 20, f = 1.7 and n = 4. The directional regularization itself has one parameter, namely the regularization weight λ. The results of the parameter optimization are presented in TableIII.• ROF -The ROF model has one parameter, namely the regularization weight λ. The results of the parameter optimization are presented in TableIII.• HDCS -The HDCS model has 4 parameter: the hybrid contrast parameter λ h , the CED contrast parameter λ c , the EED contrast parameter λ e and the scale parameter.The results of the parameters optimization are presented in TableIV. • BM4D -The original BM4D implementation has a default option that estimates the best parameters for the given image. We used this option to compute all the BM4D denoising results.

	Noise λ ROF λ directional
	2	0.001	0.010
	5	0.010	0.050
	10	0.010	0.010
	15	0.050	0.050
	20	0.010	0.050
	30	0.050	0.050
	50	0.050	0.050
		TABLE IV
		HDCS PARAMETERS.
	Noise λ h	λc	λe	scale
	2	20	40	20	0.5
	5	20	40	20	0.5
	10	40	20	20	0.5
	15	30	30	10	1.0
	20	20	40	10	1.0
	30	40	30	10	1.0
	50	40	20	10	1.0

http://vascusynth.cs.sfu.ca/Software.html

https://public.kitware.com/Wiki/TubeTK/Anisotropic Hybrid Diffusion with Continuous Switch

http://www.cs.tut.fi/ ∼ foi/GCF-BM3D

ACKNOWLEDGEMENTS

This work was funded with French Agence Nationale de la Recherche grant agreements ANR-12-MONU-0010, ANR-15-CE40-0006, ANR-18-CE45-0014 and ANR-18-CE45-0018.

APPENDIX

A. Gradient anisotropy correction


Joffrey Moiroux

Guy Boivin

Jacques Brodeur

Temperature influences host-instar selection in an aphid parasitoid: a support for the relative fitness rule

come

Temperature influences host instar selection in an aphid parasitoid: support for the relative fitness rule

Joffrey Moiroux, Guy Boivin, Jacques Brodeur

Introduction

Organisms should adopt decisions that maximise their lifetime fitness, according to the environment they face. Consequently, parents should avoid laying eggs or rearing offspring in conditions where high mortality of progeny may occur [START_REF] Eriksson | Habitat selection: Do passerines nest in association with lapwings Vanellus vanellus as defence against predators?[END_REF][START_REF] Blaustein | Oviposition habitat selection in response to risk of predation in temporary pools: mode of detection and consistency across experimental venue[END_REF]. For example, parasitoid females reduce the number of hosts parasitized when predators [START_REF] Takizawa | Effect of three species of predatory ladybirds on oviposition of aphid parasitoids[END_REF] or hyperparasitoids are present in the habitat [START_REF] Höller | Enemy-induced dispersal in a parasitic wasp[END_REF][START_REF] Petersen | Hyperparasitoid volatiles as possible foraging behaviour determinants in the aphid parasitoid Aphidius uzbekistanicus (Hymenoptera: Aphidiidae)[END_REF] or avoid laying eggs in low-quality hosts [START_REF] Vinson | The general host selection behavior of parasitoid Hymenoptera and a comparison of initial strategies utilized by larvaphagous and oophagous species[END_REF]. However under some conditions, females may adopt risk-prone behaviours and reproduce in habitats where the survival of their progeny is compromised. For example, drosophila parasitoid Leptopilina heterotoma having a shortlife expectancy tend to increase self-superparasitism following a decrease in barometric pressure, a climatic parameter associated with deteriorating foraging conditions, [START_REF] Roitberg | Life expectancy and reproduction[END_REF]. The female accepts a lower immediate fitness return, a consequence of self-superparasitism, because her prospect of fitness gain later in life is reduced.

Recently, [START_REF] Giraldeau | Risk assessment and host exploitation strategies in insect parasitoids. Behavioral ecology of insect parasitoids: from theoretical approaches to field applications[END_REF] proposed a relative fitness rule for foraging parasitoids, based not only on the quality of the environment but also on female physiological condition. This rule states that "animals that do not expect to match the average expected fitness of the population will suffer a cost because they will lose representation in the next generation, perhaps even disappear. Consequently, the chances of avoiding the genetic cost for such an organism are maximised by adopting risk-prone reproduction options." Parasitoid females should thus adopt risk-prone reproduction options when expecting a deficit in reproductive success due to their intrinsic characteristics, especially body size. To our knowledge no study has yet experimentally investigated this relative fitness rule. [START_REF] Giraldeau | Risk assessment and host exploitation strategies in insect parasitoids. Behavioral ecology of insect parasitoids: from theoretical approaches to field applications[END_REF] suggested testing this prediction using female parasitoids that are phenotypically different. For instance, smaller females with lower fecundity and reduced capacity to find new host patches are expected to be more risk-prone than larger females.

The Temperature-Size Rule (Atkinson 1994: Angilletta andDunham 2003) states that ectotherms generally grow larger at lower temperature. Accordingly, parasitoids of different size can be produced at different developmental temperatures [START_REF] Colinet | Manipulation of parasitoid size using the temperature-size rule: fitness consequences[END_REF]. Females developing at high temperature are smaller, have fewer eggs at emergence and lower capacities to find new host patches; they are thus expected to engage more risk-prone behaviour than larger females that developed at cooler temperatures. Alternatively to this hypothesis, cold-developed insects are expected to have a higher metabolic rate than warm-developed individuals when tested in a common environment because of temperature compensation following the Metabolic Cold Adaptation theory [START_REF] Wohlschlag | Metabolism of an Antarctic fish and the phenomenon of cold adaptation[END_REF][START_REF] Clarke | Seasonal acclimatization and latitudinal compensation in metabolism: do they exist?[END_REF]. Cold-developed females should therefore use their energetic resources faster, live shorter [START_REF] Le Lann | Thermal plasticity of metabolic rates linked to life history traits and foraging behaviour in a parasitic wasp[END_REF]) and thereby have fewer reproduction opportunities than females developing at higher temperature when foraging at a common temperature. This alternative hypothesis implies that cold-developed females should be more risk-prone than warm-developed females because they are time-limited.

Temperature affects reproductive success of parasitoids during both larval development and adult life. Life expectancy typically decreases as temperature increases in ectotherms (e.g., [START_REF] Loeb | On the influence of food and temperature upon the duration of life[END_REF][START_REF] Moiroux | Local adaptation of a Drosophila parasitoid: habitat-specific differences in thermal reaction norms[END_REF]. Females foraging in a hot climate have less time available for oviposition, and consequently a lower expected reproductive success than females living in a cooler environment. They should thus adopt more riskprone behaviours.

Risk-prone reproduction behaviours in parasitic wasps include oviposition in the presence of superior competitors [START_REF] Vayssade | The response of life-history traits to a new species in the community: a story of Drosophila parasitoids from the Rhône and Saône valleys[END_REF], hyperparasitoids [START_REF] Höller | Enemy-induced dispersal in a parasitic wasp[END_REF] or predators [START_REF] Takizawa | Effect of three species of predatory ladybirds on oviposition of aphid parasitoids[END_REF], self-superparasitism [START_REF] Roitberg | Life expectancy and reproduction[END_REF] or oviposition in hosts suboptimal for larval development [START_REF] Henry | Does mother really know best? Oviposition preference reduces reproductive performance in the generalist parasitoid Aphidius ervi[END_REF]). This last decision may have strong consequences on reproductive success as ovipositing in suboptimal hosts can lead to an emergence rate close to zero [START_REF] Henry | Does mother really know best? Oviposition preference reduces reproductive performance in the generalist parasitoid Aphidius ervi[END_REF]) and low quality hosts should thus be avoided in optimal conditions. In this paper, we investigated the relative role of developmental and foraging temperature on acceptance of low-quality hosts by the parasitoid Aphidius ervi Haliday (Hymenoptera: Braconidae) parasitizing the potato aphid Macrosiphum euphorbiae Thomas (Hemiptera: Aphididae). For koinobiont aphid parasitoids, second and third nymphal instars usually represent high quality hosts whereas fewer progeny emerges from eggs laid in first or fourth nymphal instars [START_REF] Colinet | Host age and fitness-related traits in a koinobiont aphid parasitoid[END_REF][START_REF] Henry | Does mother really know best? Oviposition preference reduces reproductive performance in the generalist parasitoid Aphidius ervi[END_REF].

Under laboratory conditions, A. ervi females were reared at low, medium or high temperature and next allowed to oviposit on the four host nymphal instars at one of these temperatures. We expected to observe more risk-prone behaviours, i.e. higher parasitism rate in non-optimal aphid nymphal instars, in (i) females foraging at high temperature because of reduced female longevity, in (ii) warm-developed females following the body size-related prediction, and (iii) cold-developed females following the prediction based on metabolic compensation.

Material & Methods

Insect colonies

Aphidius ervi is a solitary koinobiont parasitoid of several aphid species of economic importance. A few weeks prior to the tests, an A. ervi colony was initiated with individuals bought from BioBest Canada Company (Leamington, Canada). Parasitic wasps were reared in cages on the potato aphid Macrosiphum euphorbiae at 20±1°C, 60±10% RH, 16L: 8D photoperiod. The aphid colony was established with individuals collected from potato fields in the vicinity of Québec City, Canada, and was maintained on potato plants, Solanum tuberosum Linnaeus, cultivar "Norland" under the same conditions.

Preliminary experiment

To confirm that observation of a sting is relevant to evaluate host-instar acceptance by A. ervi, we determined if the probability of a sting resulting in the deposition of an egg was similar between host nymphal instars, as observed by [START_REF] Henry | Does mother really know best? Oviposition preference reduces reproductive performance in the generalist parasitoid Aphidius ervi[END_REF]. Five 48h-old A. ervi mated females were each allowed to attack eight aphids from each four nymphal instar. Aphids were removed after one sting and reared on a potato leaf for 48h at 20±1°C, 60±10% RH, 16L: 8D. After this period, 25 aphids per nymphal instar were randomly selected and dissected in a drop of Ringer's solution to check for parasitoid eggs. There were significantly fewer aphids that contained eggs than aphids that were stung (z = 3.114, p < 0.001), with 65 to 70% of attacked aphids being parasitized. However, this pattern was similar across host instar (GLM, z = 0.608, p = 0.517).

Observation of a sting can therefore be used to assess host-instar acceptance by A. ervi.

Experimental design

To measure the influence of developmental and foraging temperature on host selection behaviour in A. ervi, we first reared parasitoids at three different temperatures: 12, 20 and 28±1°C. These temperatures were selected because 20°C is the optimal temperature for rearing while 12 and 28°C are close to the extremes at which A. ervi remains active when foraging [START_REF] Gilchrist | A quantitative genetic analysis of thermal sensitivity in the locomotor performance curve of Aphidus ervi[END_REF]. Fifty parasitoid couples were each placed in a Petri dish at 20±1°C for five hours with 30 N2 aphids placed on a potato leaf. After this period, 10 presumably parasitized aphids per parasitoid couple were reared at 12, 20 and 28±1°C on a potato plant. Once parasitized aphids turned into mummies, they were placed individually in gelatine capsules until parasitoid emergence. Emerged female parasitoids were then kept for 24 hours in a Petri dish containing moistened cotton and honey along with a male from the same treatment at the same temperature they experienced during development. This period allowed females to mate and mature eggs (Le Lann et al.

2010).

The next day, the host selection behaviour of A. ervi females was characterized for each of the four aphid nymphal instars (N1 to N4) in a climate room (60±10% RH). Seven treatments were established to examine the relative role of developmental and foraging temperatures on host selection: (1-3) females foraging at the same temperature as the one they experienced during their development (i.e., 12, 20 and 28°C; named 12-12, 20-20 and 28-28 treatments, respectively); (4-5) females reared at 20°C and foraging at 12 or 28°C (20-12 and 20-28 treatments, respectively); (6-7) females reared at 12°C and 28°C and foraging at 20°C (12-20 and 28-20 treatments, respectively).

One hour before the beginning of the observation, a potato leaf (ø 4cm) was taped on a piece of paper and 60 apterous aphids were placed on it, 15 from each nymphal instar.

Concurrently, parasitoid females were acclimatized to the experimental temperature for one hour. Females were next introduced in the experimental arena (ø 26cm, h 8.5 cm) and their behaviour was monitored using a camera (Panasonic WV-BD400) and next analyzed with The Observer XT software (Noldus, Wageningen, The Netherlands).

Twenty-five females per treatment were observed and the frequency of the following behaviours was determined: successful oviposition (sting), failed oviposition (unsuccessful oviposition attempt due to aphid behavioural defence), host rejection (parasitoid contact with a motionless aphid not followed by an attack) and misplaced attack (sting attempt on a section of the potato leaf with no aphid). This latter behaviour is rarely taken into account in similar studies but was frequently observed in some treatments and could indicate a neural disruption induced by temperature. Handling time was measured and included duration of host antennation, oviposition and time to recover from behavioural defence. We did not replace parasitized aphids because superparasitism was rare under our experimental conditions (<5%, J. Moiroux, unpublished data). The numbers of successful ovipositions, failed ovipositions and host rejections were pooled to calculate the number of encounters per host nymphal instar. An attack rate per host nymphal instar was also calculated, i.e. the sum of successful and failed ovipositions on a given nymphal instar divided by the number of encounters.

The observation stopped when the parasitoid remained motionless for more than 15 minutes or left the patch for more than one minute at 28°C, 2 minutes at 20°C and 6 minutes at 12°C, these durations being proportional to the relationship between temperature and walking speed for A. ervi females [START_REF] Gilchrist | A quantitative genetic analysis of thermal sensitivity in the locomotor performance curve of Aphidus ervi[END_REF]. At the end of the experiment, parasitoid females were frozen and their tibia length was measured under a microscope with the µScope Microscopy Software (PixeLINK®, Ottawa, Canada). Hosts from a given treatment were sorted by nymphal instar, placed on a potato plant and reared at 20±1°C, 60±10% RH, 16L: 8D photoperiod until parasitoid emergence.

Mummification rate, emergence rate and tibia length of A. ervi progeny were determined to assess the suitability of each aphid nymphal instar for parasitoid larval development.

We only considered tibia length of female progeny because no males emerged from N4 aphids and a sexual dimorphism in body size exists in A. ervi [START_REF] Sequeira | Seasonal variation in body size and offspring sex ratio in field populations of the parasitoid wasp, Aphidius ervi (Hymenoptera: Aphidiidae)[END_REF].

Longevity

Females reared at 12, 20 and 28°C and not used in the previous experiment were placed at emergence at 20°C, 60±10% RH, 16L: 8D in 12 cm cubic cages containing moistened cotton and honey to measure their longevity. Mortality was recorded twice a day, in the morning and evening.

Statistical analyses

The effect of both developmental temperature and foraging temperature on body size of parasitoid females was analyzed with a two-way ANOVA, data being normally distributed.

We compared the number of successful and failed ovipositions and the number of host rejections using Generalized Linear Models (GLM) with developmental temperature, foraging temperature and host nymphal instar as fixed factors and the number of encounters as a covariate, considering a Poisson distribution of data. Attack rate was compared using a Generalized Linear Model (GLM) with developmental temperature, foraging temperature and host nymphal instar as fixed factors, considering a binomial distribution of data. The number of misplaced attacks was analysed using a GLM with developmental and foraging temperature as fixed factors. We also performed GLM with developmental temperature, foraging temperature and host nymphal instar as factors to test for differences in handling time, mummification rate, emergence rate and progeny tibia length. The survival package provided by R software was used to test for differences in longevity between females that developed at different temperatures, using a Weibull distribution. All statistical analyses were carried out using R software version 2.14.1 (R Development Core Team, 2011).

Results

Developmental temperature significantly influenced A. ervi female size (F= 7.14, p < 0.001), tibia length decreasing with increasing temperature (12°C = 0.76 ± 0.06 mm, 20°C = 0.67 ± 0.05 mm, 28°C = 0.56 ± 0.05). In contrast, no difference in parasitoid size was associated to foraging temperature (F= 0.49, p = 0.806).

The attack rate was significantly influenced by host nymphal instar (z = 4.118, p = 0.012), developmental temperature (z =5.907, p < 0.001), foraging temperature (z = 6.252, p < 0.001), the interaction between nymphal instar and developmental temperature (z = 3.816, p = 0.031) and the interaction between nymphal instar and foraging temperature (z = 3.057, p = 0.043). Second and third nymphal instars were more frequently attacked than first and fourth nymphal instars. Females that developed at 28°C attacked first and fourth nymphal instars more frequently than females that developed at 12 and 20°C, and females foraging at 28°C attacked N4 aphids more frequently than females foraging at lower temperatures.

Host stage (z = 3.383, p < 0.001), developmental temperature (z = 2.388, p = 0.017) and the interaction between developmental temperature, foraging temperature and host nymphal instar (z = -3.47, p < 0.001) significantly influenced the number of successful ovipositions by A. ervi. Second and third nymphal instars were more frequently parasitized than other instars. First nymphal instars were more frequently parasitized by females that developed at 28°C than females that developed at 12°C (p = 0.018) or 20°C (p < 0.001), with no difference between these two treatments (p = 0.617). The numbers of second and third nymphal instars parasitized by A. ervi were the highest in the 20-20 and 20-28 treatments and the lowest in the 28-28 treatment. Fourth nymphal instars were the most frequently parasitized in the 20-28 treatment, second in the 12-12 treatment and no difference was observed between the other treatments (Figure 1A).

The number of failed ovipositions varied with host age (z = 2.573, p = 0.010) and the interaction between developmental and foraging temperature (z = 2.457, p = 0.014).

Parasitoids failed more frequently to parasitize N4 aphids than other host ages and females from the 28-28 treatment were less successful than other females (Figure 1C), especially when attacking N4 aphids.

Host nymphal instar (z = 5.758, p <0.001), developmental temperature (z = 6.559, p < 0.001) and the interaction between developmental temperature and host nymphal instar (z = -2.803, p = 0.005) significantly influenced the number of host rejections. First nymphal instars were more frequently rejected than other host nymphal instars in all treatments.

Rejections of N1 were most frequent in females that developed at 20°C and less frequent in females that developed at 28°C (Figure 1B).

The number of misplaced attacks was significantly influenced by developmental temperature (z = 3.917, p < 0.001), foraging temperature (z = 3. 514, p < 0.001) and the interaction between these two factors (z = 4.712, p <0.001). It was the highest in the 12-12 treatment and higher for females that developed or foraged at 12°C than for females from other treatments (Number of misplaced attacks, mean ± se; 12-12: 9.9 ± 0.9; 20-20:

1.0 ± 0.4; 28-28: 1.0 ± 0.5; 20-12: 5.7 ± 0.6; 20-28: 0.7 ± 0.3; 12-20: 4.4 ± 0.7; 28-20: 1.7 ± 0.4).

Host handling time was influenced by host nymphal instar (z = -2.751, p = 0.027), and the interaction between nymphal instar, developmental and foraging temperature (z = 3.989, p < 0.001). Parasitoids took more time to attack N4 aphids than other instars, N1 aphids being attacked the most quickly. Females from the 20-28 treatment handled hosts faster than other females, while handling time was the longest for females from the 28-28 treatment (Figure 2). Handling time decreased with developmental temperature for females foraging at 20°C.

Development

Mummification rate (z = 3.761, p < 0.001) and emergence rate (z = 3. 012, p = 0.003) of A. ervi were only influenced by host age (Table 3), being highest in N2 aphids, intermediate in N3 and lowest in N1. Size of A. ervi progeny was influenced by host nymphal stage (F = 4.117, p = 0.012). Parasitoids laid in N3 aphids were bigger than parasitoids laid in other nymphal instars, with parasitoids laid in N1 aphids being the smallest (Table 1).

Longevity

Developmental temperature significantly influenced longevity of parasitoid females (Survival analysis, χ2 = 48.07, p = 0.013). When transferred to 20°C at emergence, wasps that developed at 12°C (longevity = 8.76•± 1.48 days, mean•± SE) lived significantly shorter than wasps that developed at 20°C (15.43•± 2.17 days) and 28°C (12.86•± 1.93 days), with no difference being observed between the last two treatments (p = 0.247).

Discussion

Temperature affects body size, initial egg load and basal metabolic rate during A. ervi larval development, as well as their life expectancy during adult life. Following the relative fitness rule proposed by [START_REF] Giraldeau | Risk assessment and host exploitation strategies in insect parasitoids. Behavioral ecology of insect parasitoids: from theoretical approaches to field applications[END_REF], we predicted that (1) warm-developed parasitoid females should adopt more risk-prone behaviours than females that developed at lower temperature because of their smaller body size and egg load; (2) cold-developed females should adopt more risk-prone behaviours than females developed at warmer temperature when foraging in a common environment because of metabolic compensation; (3) females foraging at high temperature should adopt more risk-prone behaviours than females foraging at lower temperature because of reduced longevity. The first and third predictions are supported by our results as A. ervi (1) females that developed at high temperature attacked low-quality hosts more frequently than other females and (2) females foraging at high temperature attacked low-quality hosts more frequently than those foraging at lower temperature. The second prediction was not supported, as females that developed at cold temperature attacked similar proportions of less suitable hosts than females that developed at intermediate temperature.

Optimal thermal conditions

When A. ervi females developed and foraged at 20°C, their optimal thermal conditions, they preferentially parasitized second and third nymphal instar aphids and avoided first and fourth nymphal instars. [START_REF] He | Host age preference in Aphidius ervi (Hymenoptera: Aphidiidae)[END_REF] also observed a preference for intermediary nymphal instars for A. ervi parasitizing the pea aphid Acyrthosiphon pisum.

We observed that mummification rate, emergence rate and body size of progeny were the highest when A. ervi developed in N2 and N3 hosts. Host selection by A. ervi females is thus optimal at 20°C. Avoidance of first and fourth nymphal instar may result from a combination of factors including lower nutritional quality of small hosts [START_REF] Vinson | Host suitability for insect parasitoids[END_REF], increased resistance in larger/older hosts [START_REF] Li | Pea aphid clonal resistance to the endophagous parasitoid Aphidius ervi[END_REF] and more efficient behavioral defence by fourth instar aphids that lower their profitability to parasitoids [START_REF] Barrette | Testing competing measures of profitability for mobile resources[END_REF].

Developmental temperature

We observed that female parasitoids developing at high temperature (28-28 and 28-20 treatments) and having a lower initial egg load (Moiroux et al. submitted), increased their attack rate on first and fourth nymphal instars, despite an emergence rate ranging from 2 to 4% in these hosts. Females thus took the risk of ovipositing in hosts where their progeny faced a high mortality rate. The higher number of failed ovipositions in N2 to N4 instars observed for 28-28 females may partly explain higher acceptance of N1 hosts. However, 28-20 females succeeded as much as females that developed at lower temperatures to oviposit in N2 to N4 aphids but still attacked N1 aphids. This suggests that females attacked first nymphal instar not because they were unable to parasitize other instars but rather to increase the number of oviposition opportunities. This result is in contradiction with the prediction that egg-limited parasitoids, i.e. warm-developed females in our study, should accept only high quality hosts for oviposition [START_REF] Mangel | Evolution of host selection in parasitoids: does the state of the parasitoid matter[END_REF][START_REF] Heimpel | Egg limitation, host quality, and dynamic behavior by a parasitoid in the field[END_REF]. Time-limitation, which is expected to result in acceptance of low-quality hosts [START_REF] Mangel | Evolution of host selection in parasitoids: does the state of the parasitoid matter[END_REF], is also unlikely to explain our results since 28-20 females lived as long as 20-20 females, and longer than 12-20 females. The relative fitness rule [START_REF] Giraldeau | Risk assessment and host exploitation strategies in insect parasitoids. Behavioral ecology of insect parasitoids: from theoretical approaches to field applications[END_REF] however helps to explain our results.

Warm-developed females were smaller than females from other treatments, emerged with a lower initial egg load and less energetic resources required for egg maturation [START_REF] Giron | Lipogenesis in an adult parasitic wasp[END_REF] and foraging activity (Nation 2008). According to the fitness rule, these small females may expect a lower reproductive success than the overall population and should therefore adopt risk-prone behaviours, i.e. oviposition in less suitable hosts. A small female is expected to have a lower fitness than a large female because of lower fecundity and smaller capacity to find and attack hosts. To increase their relative fitness, an optimal decision woud be to parasitize a maximum number of hosts as rapidly as possible, relying on occurrence of a sudden high mortality event, such as a heat wave or a storm.

The lower attack rate of low-quality hosts by females from the 12-20 treatment in comparison to females that developed at higher temperature and foraged at 20°C does not support the alternative prediction that temperature compensation in cold-developed organisms [START_REF] Wohlschlag | Metabolism of an Antarctic fish and the phenomenon of cold adaptation[END_REF][START_REF] Clarke | Seasonal acclimatization and latitudinal compensation in metabolism: do they exist?[END_REF]) should increase risk-prone behaviours. The shorter longevity and handling time we observed for females that developed at low temperature both suggest that metabolic compensation occurs in A. ervi, as observed for the closely related species Aphidius rhopalosiphi [START_REF] Le Lann | Thermal plasticity of metabolic rates linked to life history traits and foraging behaviour in a parasitic wasp[END_REF]. It is thus unlikely that the pattern observed for 12-20 females resulted from the absence of temperature compensation. If parasitoids assess their expected reproductive success through egg load and life expectancy, the higher initial egg load in females that developed at low temperature may compensate for the lower life expectancy when foraging at intermediate temperature and may thus influence risk-sensitive decision-making.

Foraging temperature

Developmental temperature clearly influenced oviposition decisions of parasitoids foraging at a common temperature. Foraging temperature also strongly influenced host selection by A. ervi females. Parasitoids foraging at high temperature (20-28 and 28-28 treatments) attacked low-quality hosts more frequently than adult females experiencing lower temperature. They attacked N4 aphids at a similar rate and more frequently than in other treatments but 20-28 females were more successful than 28-28 females as they were fast enough to avoid kicks from N4 aphids. Although both foraging at 28C, females from these two treatments experienced different profitability [START_REF] Barrette | Testing competing measures of profitability for mobile resources[END_REF]) because handling time and oviposition success varied with developmental temperature. These results support the hypothesis that foraging females associate high temperature to short life expectancy and thus low reproductive success, resulting in risk-prone behaviours. [START_REF] Roitberg | Life expectancy and reproduction[END_REF] observed a similar result in Leptopilina heterotoma females which superparasitized more frequently when exposed to a probability of decreased life expectancy, simulated by a rapid drop in atmospheric pressure as observed before a storm.

Aphidius ervi females from the 12-12 treatment rejected hosts very rarely compared to females from all other treatments. Furthermore, the number of misplaced attacks was the highest in the 12-12 treatment. This suggests that an alteration of the host selection behaviour, including host-instar discrimination, occurred at low temperature, as observed in cold stored Anaphes victus females (van Baaren, [START_REF] Van Baaren | Effect of low temperature exposure on oviposition behaviour and patch exploitation strategy in parasitic wasps[END_REF]. Effects of cold temperature on nervous system would explain the acceptance of all aphid nymphal instars by 12-12 females although we cannot reject the hypothesis of a behavioural change with low temperature. Developmental and foraging temperatures may be both involved in this alteration, as we observed an increased number of misplaced attacks in 12-20 and 20-12 treatments.

PerformanceOur experimental design allowed us to test the classical physiological hypothesis that performance is maximal at the temperature at which an organism is acclimated. Our results concur with other studies in which no beneficial effect was associated to acclimation (e.g., [START_REF] Leroi | Temperature acclimation and competitive fitness: an experimental test of the beneficial acclimation hypothesis[END_REF][START_REF] Gibert | Locomotor performance of Drosophila melanogaster: interactions among developmental and adult temperature, age, and geography[END_REF], as females from the 12-12 and 28-28 treatments parasitized fewer hosts than females from other treatments. Our results however support the optimal developmental temperature hypothesis [START_REF] Cohet | Control of the adult reproductive potential by preimaginal thermal conditions: a study in Drosophila melanogaster[END_REF] that states that ectotherms that developed at intermediate temperatures perform well at a wider range of temperatures than individuals developed at extreme temperatures. Aphidius ervi females reared at 20°C indeed performed better than females from other treatments, independently of foraging temperature, as observed by [START_REF] Gibert | Locomotor performance of Drosophila melanogaster: interactions among developmental and adult temperature, age, and geography[END_REF] in Drosophila melanogaster.

Conclusion

An increasing number of studies have identified potential consequences of global warming on trophic relationships and the structure and stability of communities.

Temporal and spatial species associations between trophic levels may be sensitive to the expected increase in mean temperature and frequency of extreme climatic events [START_REF] Parmesan | Ecological and evolutionary responses to recent climate change[END_REF]. Functional and numerical responses have also been shown to be temperature-dependent (e.g., [START_REF] Sentis | Using functional response modeling to investigate the effect of temperature on predator feeding rate and energetic efficiency[END_REF] and other behavioural components are should be taken into account when investigating consequences of climate change on interacting species. We present evidence that developmental and foraging temperatures are both involved in host selection and risk-sensitive decision-making in parasitoids. For instance, an increase in developmental and/or foraging temperature resulted in acceptance of less profitable hosts, which could have detrimental consequences on parasitoid populations in natural and managed ecosystems. Our study underlines the necessity to further investigate the influence of temperature on behaviour to understand how organisms and communities will cope with global warming. 4. We observed that parasitoids foraging at high temperature attacked low-quality hosts (first and fourth nymphal instars) more frequently than females from other treatments and that warm-developed females attacked first nymphal instars more frequently than females from other treatments. These results support our hypothesis that a decrease in the expected parasitoid reproductive success, indicated by a lower life expectancy or initial egg load, resulted in risk-prone behaviours. Cold-developed females did not attack nonoptimal hosts more frequently than females. This result does not support the hypothesis that higher metabolic rate results in risk-prone behaviours.

5. To our knowledge, this is the first experimental evidence that temperature influences host stage selection and risk-sensitive making-decision in parasitoids and the first study to support the relative fitness rule.

Keywords: Optimal foraging, climate change, host selection, Aphidius ervi,

Figure 1 .

 1 Figure 1. Mean (± s.e) number of (A) successful ovipositions (B), failed ovipositions and (C) host rejections in each of the four nymphal instar of Macrosiphum euphorbiae by Aphidius ervi females exposed to seven thermal treatments.

Figure 2 .

 2 Figure 2. Mean (± s.e) handling time of Aphidius ervi females attacking Macrosiphum euphorbiae aphid hosts of four different nymphal instars in seven thermal treatments.

Table 1 .

 1 Mummification rate, emergence rate and tibia length of female progeny measured on parasitoids from seven thermal treatments where Aphidius ervi females were allowed to parasitize 60 Macrosiphum euphorbiae aphids of four different instars. Mean ± s.e.

			Host stage
		N1	N2	N3	N4
	Mummification rate	0.05 ± 0.03 a 0.31 ± 0.08 b 0.23 ± 0.06 b 0.08 ± 0.04 a
	Emergence rate	0.02 ± 0.01 a 0.27 ± 0.05 b 0.16 ± 0.04 c 0.04 ± 0.02 a
	Tibia length of female progeny (mm)	0.51 ± 0.03 a 0.64 ± 0.04 b 0.69 ± 0.05 c 0.56 ± 0.03 d

Acknowledgments

We are grateful to Josée Doyon for technical assistance during this study. This research was part of the project "Impact of climate change on synchronism between pests and their natural enemies" supported by the Consortium on Regional Climatology and Adaptation to Climate Change, OURANOS (Québec, Canada), and the Fonds vert of the Québec Ministry of Agriculture, Fisheries and Food.