Maroua Mehri
email: maroua.mehri@univ-lr.fr

Nabil Sliti
email: selliti.nabil@gmail.com

Pierre Héroux
email: pierre.heroux@univ-rouen.fr

Petra Gomez-Krämer

Najoua Essoukri

Ben Amara

Rémy Mullot
email: remy.mullot@univ-lr.fr

Use of SLIC superpixels for ancient document image enhancement and segmentation

Keywords: Ancient document images, enhancement, segmentation, SLIC superpixels, Gabor features, multiscale approach

Designing reliable and fast segmentation algorithms of ancient documents has been a topic of major interest for many libraries and the prime issue of research in the document analysis community. Thus, we propose in this article a fast ancient document enhancement and segmentation algorithm based on using Simple Linear Iterative Clustering (SLIC) superpixels and Gabor descriptors in a multi-scale approach. Firstly, in order to obtain enhanced backgrounds of noisy ancient documents, a novel foreground/background segmentation algorithm based on SLIC superpixels, is introduced. Once, the SLIC technique is carried out, the background and foreground superpixels are classified. Then, an enhanced and non-noisy background is achieved after processing the background superpixels. Subsequently, Gabor descriptors are only extracted from the selected foreground superpixels of the enhanced gray-level ancient book document images by adopting a multi-scale approach. Finally, for ancient document image segmentation, a foreground superpixel clustering task is performed by partitioning Gabor-based feature sets into compact and well-separated clusters in the feature space. The proposed algorithm does not assume any a priori information regarding document image content and structure and provides interesting results on a large corpus of ancient documents. Qualitative and numerical experiments are given to demonstrate the enhancement and segmentation quality.

INTRODUCTION

Since the early 1990s, libraries and museums have conducted large digitization campaigns with cultural heritage documents and scientific resources for ensuring a restoration and lasting preservation of historical collections and promoting worldwide accessibility to cultural patrimony which requires to be protected from further deterioration and damages caused by repetitive handling [START_REF] Baird | Digital libraries and document image analysis[END_REF]. Due to the huge amount of numeric high-quality reproductions induced by the rapid growth of digital libraries worldwide, many challenges and open issues have been raised, such as providing a reliable document interpretation system and developing an efficient content-based image retrieval tool. Thus, a raising interest to ancient document image analysis has been generated since it helps to reach the objective of ensuring a relevant indexing and retrieval of digitized resources and offering a structured access to large sets of cultural heritage documents [START_REF] Coustaty | Historical document analysis: A review of French projects and open issues[END_REF]. The ancient document image analysis starts by segmenting a document in order to find and classify homogeneous regions or zones, such as graphic and textual regions [START_REF] Okun | A survey of texture-based methods for document layout analysis[END_REF]. Finding graphic regions can be used to segment and analyze the graphical part in historical heritage, such as the drop caps (i.e., the first letters at the beginning of a paragraph) [START_REF] Uttama | Segmentation and retrieval of ancient graphic documents[END_REF] while determining text zones can be used as a pre-processing stage for character recognition [START_REF] Diem | Recognition of degraded handwritten characters using local features[END_REF], text line extraction [START_REF] Likforman-Sulem | Text line segmentation of historical documents: A survey[END_REF] and handwriting recognition [START_REF] Fischer | A combined system for text line extraction and handwriting recognition in historical documents[END_REF], etc. In this work, we are interested in ancient document image analysis. Our goal consists of identifying homogenous regions or similar groups of pixels sharing some visual characteristics by labeling and grouping pixels from ancient document images.

Recently, a family of pixel-based algorithms for document image content and region extraction in document images has been proposed [START_REF] Baird | Document content inventory and retrieval[END_REF][START_REF] An | Iterated document content classification[END_REF][START_REF] An | Document segmentation using pixel-accurate ground truth[END_REF]. The document image analysis issue had been considered as texture segmentation and classification [START_REF] Okun | A survey of texture-based methods for document layout analysis[END_REF]. Kise precised that the most relevant methods used to analyze pages with overlapping or unconstrained layouts are based on signal properties of page components by investigating texturebased features and techniques [START_REF] Kise | Page segmentation techniques in document analysis[END_REF]. Thus, in this article a texture-based approach is proposed in order to assist the document image analysis without hypothesis on the document structure, neither on the document model nor the typographical parameters. To our knowledge, there are a very limited number of comparative studies of texture-based approaches in the field of ancient document characterization and segmentation. For instance, Mehri et al. compared some of the well-known texture-based approaches: the auto-correlation function [START_REF] Petrou | Image processing: Dealing with texture[END_REF], Gray Level Co-occurrence Matrix (GLCM) [START_REF] Haralick | Textural features for image classification[END_REF], and Gabor filters [START_REF] Gabor | Theory of communication. Part 1: The analysis of information[END_REF], used for segmentation of digitized simplified ancient document images [START_REF] Mehri | Texture feature evaluation for segmentation of historical document images[END_REF]. Qualitative and quantitative results were obtained on simplified ancient document images. They concluded that Gabor features performed better than auto-correlation and GLCM ones for font segmentation on the one hand, and for distinguishing textual regions from graphical ones on the other hand. Therefore, Gabor features are investigated in this article for ancient document image analysis.

Nevertheless, the literature shows that many different problems may arise related to the particularities of ancient documents such as a large variability of the page layout: noise and degradation caused by copying, scanning and aging (staining, mold or moisture, and the faded out ink and uneven lighting due to the folded and corrugated parchment or papyrus, etc.), page skew, complicated layout, random alignment, specific fonts, presence of embellishments, irregular spacings between the characters, words, lines, paragraphs and margins, overlapping object boundaries, varying text column widths, interspersed graphics, text in different orientations, presence of black borders, and superimposition of information layers (stamps, handwritten notes, noise, back-tofront interference or ink that was bleeding through, etc.). These characteristics complicate the characterization or segmentation of ancient documents and make the processing of this kind of document a difficult task [START_REF] Agam | Degraded document image enhancement[END_REF][START_REF] Coustaty | Historical document analysis: A review of French projects and open issues[END_REF]. Thus, a pre-processing phase is necessary to remove background noise (e.g. superfluous information appearing from the verso, etc.). Indeed, in this work a novel foreground-background segmentation algorithm is introduced before the Gabor feature extraction and analysis task to obtain enhanced backgrounds of noisy ancient documents. The proposed algorithm is based on Simple Linear Iterative Clustering (SLIC) superpixels. Therefore, a fast automatic ancient document enhancement and segmentation algorithm based on using SLIC superpixels and Gabor descriptors in a multi-scale approach is proposed in this article.

The remainder of this article is organized as follows: in Section 2, the superpixel technique used for document image enhancement is described briefly. Section 3 details different Gabor-based approaches proposed in the literature and presents the investigated Gabor features. Section 4 presents the proposed ancient document enhancement and segmentation algorithm. In Section 5, we outline the experimental protocol by describing the experimental corpus and ground-truth. In addition, we discuss the pertinence of experiments by computing a classification metric for an evaluation of accuracy. Qualitative results are also given to demonstrate the enhancement and segmentation quality. Our conclusions and future work are presented in Section 6.

SUPERPIXEL TECHNIQUE

The superpixel technique has recently been used and has increasingly been applied in many computer vision applications, such as Gaussian noise estimation [START_REF] Wu | Gaussian noise estimation with superpixel classification in digital images[END_REF], class segmentation and object localization [START_REF] Fulkerson | Class segmentation and object localization with superpixel neighborhoods[END_REF], etc. Instead of using a rigid structure of pixel grid for feature extraction and processing at the pixel level for segmentation, localization, and classification issues, the superpixel technique has been used as the basic unit and become a consistent alternative. The superpixel approach has the advantage to be faster, more memory efficient, and more interesting to compute image features on each superpixel center than on each image pixel [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. It consists of grouping pixels sharing similar characteristics or properties (e.g. texture cues, contour, color, etc.) into a significant polygon-shaped region. Indeed, the superpixel approach produces an over-segmented image representing a compact content map [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF].

Recently, the superpixel approach has gained great attention of many researchers in document image analysis fields. For instance, Cohen et al. separated drawings from background and noise of ancient documents by using spatial and color features which were extracted from superpixels [START_REF] Cohen | Robust text and drawing segmentation algorithm for historical documents[END_REF]. Mehri et al. used an entropy rate superpixel segmentation algorithm based on a novel objective function [START_REF] Liu | Entropy rate superpixel segmentation[END_REF] to obtain enhanced backgrounds of noisy ancient documents [START_REF] Mehri | Robustness assessment of texture features for the segmentation of ancient documents[END_REF]. By setting the number of superpixels equal to 5% of image pixels, the entropy rate superpixel segmentation algorithm is carried out on ancient document image. Then, the background and foreground superpixels are classified. Subsequently, an enhanced and non-noisy background is achieved after processing the background superpixels, i.e. the background superpixels and the pixels belonging to them have a value of white pixel (i.e. a 255 gray-level value).

Many superpixel algorithms have been proposed in the literature. Achanta et al. classsified the existing superpixel methods into three classes: the graph-based, gradient-ascent-based, and SLIC superpixels [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. The graph-based approaches consider each image pixel as a node in a graph. Then, the superpixels are generated by minimizing a cost function based on a the obtained graph. The gradient-ascent-based approaches process by iteratively refine the classification of samples on homogeneous and compact clusters until a convergence criterion is achieved to form superpixels. Finally, the SLIC approach is based on grouping pixels based on their color similarity and proximity in the image plane, by considering the five-dimensional space [l, a, b, x, y], where {l, a, b} represent the three pixel coordinates in the CIELAB color space, and {x, y} are the pixel spatial coordinates.

The SLIC technique is based on an adaptive form of k-means clustering approach to well produce superpixels. By initiating the desired number of superpixels k s , the SLIC algorithm defines the k s initial cluster centers

C i = [l i , a i , b i , x i , y i ] T
, where 0 ≤ i ≤ k s , which form a structure of regular equally-sized superpixel grid. In order to generate superpixels with similar sizes, the grid interval S is set to N ks , where N represents the number of pixels. Then, each pixel is assigned to the nearest cluster center. Finally, an update task is performed to control the cluster center according to the mean [l, a, b, x, y] T vector of all the pixels belonging to the cluster. Achanta et al. suggested to repeat 10 times the assignment and update steps of the SLIC approach to ensure good segmentation results [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. They stated the SLIC technique as the best superpixel method comparing to five state-of-the-art superpixel methods in terms of the adherence to boundaries, computational and memory efficiency, segmentation performance, ease of use, parameter tuning, ability to control superpixel compactness, regularity and number, etc.

GABOR FEATURES

The investigated texture features tested in our work are Gabor descriptors which are extracted using the multichannel Gabor filtering technique. The original Gabor elementary functions were firstly proposed by Gabor [START_REF] Gabor | Theory of communication. Part 1: The analysis of information[END_REF]. The multi-channel Gabor filtering was inspired by the multi-channel filtering theory first proposed by Campbell and Robson for processing of visual information of the human visual system [START_REF] Campbell | Application of Fourier analysis to the visibility of gratings[END_REF]. Daugman modeled the processing of visual information of the human visual system by the two-dimensional multi-channel Gabor functions, which are local spatial bandpass filters [START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters[END_REF]. The main idea of the multi-channel filtering technique is to exploit the differences in dominant sizes and orientations of different textures by decomposing the original image into several filtered images with limited spectral information. The two-dimensional Gabor functions have the advantage to have conjoint resolution information in both the two-dimensional spatial and Fourier domains. The filtered images are proceeded by tuning the analyzed image to combinations of frequency and orientation in a narrow range, which are referred to channels and interpreted as band-pass filters. By applying a bank of Gabor filters, the specified channels cover the spatial-frequency domain. Ursani et al. presented an empirical comparison between texture features based on the discrete Fourier transform and Gabor filters for texture recognition and retrieval [START_REF] Ursani | Texture features based on Fourier transform and Gabor filters: An empirical comparison[END_REF]. They proved that analyzing Gabor features in image datasets containing noisy and rotated variants of texture performed better than analyzing Fourier descriptors for texture recognition and retrieval. Hence, Gabor filters have been shown to have good performance, due to its optimal localization properties to capture information in both spatial and frequency domains from the analyzed images, as opposed to Fourier transform.

A two-dimensional Gabor filter is a linear selective band-pass filter, dependent on two parameters: spatial frequency f g and orientation θ g , which characterize the specified channel. It consists of a Gaussian kernel function modulated by a sinusoidal plane wave. The spatial frequency f determines the distance from the Gaussian centers to the origin while the orientation θ g specifies the angle from the horizontal axis, i.e. α-axis to the Gaussian centers. The multi-channel Gabor filtering approach is inherently multi-resolutional, which is a close relative of the wavelet transform [START_REF] Jain | Unsupervised texture segmentation using Gabor filters[END_REF].

The Gabor transform of an image I(x, y) is:

I G (fg ,θg ) (x, y) = α∈Ω β∈Ω I(x + α, y + β) * G (fg,θg) (α, β) (1)
where f g , θ g and σ are the spatial frequency, orientation, and space constant of the Gabor envelope.

Texture features generated by Gabor filters have been increasingly considered and applied to document image analysis. During the last two decades, Gabor texture analysis approaches have been proposed for biometric identification based on handwriting [START_REF] Zhu | Biometric personal identification based on handwriting[END_REF][START_REF] Feddaoui | Personal identification based on texture analysis of Arabic handwriting text[END_REF][START_REF] Chen | Gabor features for offline Arabic handwriting recognition[END_REF], writer identification [START_REF] Bensefia | A writer identification and verification system[END_REF], handwritten word recognition [START_REF] Buse | A structural and relational approach to handwritten word recognition[END_REF], character recognition [START_REF] Wang | Gabor filters-based feature extraction for character recognition[END_REF], font recognition [START_REF] Zhu | Font recognition based on global texture analysis[END_REF], script identification [START_REF] Tan | Rotation invariant texture features and their use in automatic script identification[END_REF][START_REF] Joshi | Script identification from Indian documents[END_REF], signature recognition [START_REF] Kekre | Gabor filter based feature vector for dynamic signature recognition[END_REF], palm print recognition [START_REF] Mu | Mean and standard deviation as features for palmprint recognition based on Gabor filters[END_REF], etc. Zhu et al. proposed a texture-analysis-based algorithm for automatic font recognition by extracting Gabor features [START_REF] Zhu | Font recognition based on global texture analysis[END_REF]. The authors showed a 99, 1% of mean recognition rate. Ma and Doermann proposed a Gabor filter based multi-class classifier in order to identify scripts, font-faces and font-styles [START_REF] Ma | Gabor filter based multi-class classifier for scanned document images[END_REF].

Nevertheless, numerous approaches have been sought for text segmentation and extraction from digital document images using Gabor descriptors [START_REF] Jain | Text segmentation using Gabor filters for automatic document processing[END_REF][START_REF] Chan | Text analysis using local energy[END_REF][START_REF] Raju | Gabor filters for document analysis in Indian bilingual documents[END_REF][START_REF] Raju | Text localization and extraction from complex color images[END_REF]. Several studies have been conducted in the literature for page layout analysis using the multi-channel Gabor filters [START_REF] Randen | Segmentation of text/image documents using texture approaches[END_REF][START_REF] Ma | Gabor filter based multi-class classifier for scanned document images[END_REF][START_REF] Pati | Gabor filters for document analysis in Indian bilingual documents[END_REF], while few ones have explored Gabor filters for ancient document image segmentation. For instance, Ribeiro et al. proposed an optical character recognition system for ancient document recognition by applying fuzzy methods on aligned oriented features extracted using Gabor filters in the training step [START_REF] Ribeiro | Ancient document recognition using fuzzy methods[END_REF]. Vieux and Domenger proposed a pixel-based classification approach to separate text from other classes (e.g. illustrations and background) by using a bank of Gabor filters at five scales (1, √ 2, 2, 2 √ 2 and 4) and six orientations (k π 6 , k ∈ {0, . . . , 5}) [START_REF] Vieux | Hierarchical clustering model for pixel-based classification of document images[END_REF]. They evaluated their approach on a public dataset containing magazines and technical journals. They found 86%, 82.7%, and 53.7% of F-measure for segmenting background, text, and graphic pixels, respectively. Jain et al. showed the effectiveness of applying a multi-channel Gabor filtering-based texture segmentation approach for the segmentation and classification of document images [START_REF] Jain | Page segmentation using texture analysis[END_REF]. They chose the five following spatial frequencies: 4

√ 2, 8 √ 2, 16 √ 2, 32 √ 2 and 64 √ 2.
Charrada and Ben Amara extracted nets from ancient Arab periodicals by exploring Gabor filters [START_REF] Charrada | Texture approach for nets extraction application to old Arab newspapers images structuring[END_REF]. Zhong and Cheriet used the dimensionally reduced multi-channel Gabor filters for text block identification on image patches from ancient documents [START_REF] Zhong | Image patches analysis for text block identification[END_REF]. The authors extracted [START_REF] Feddaoui | Personal identification based on texture analysis of Arabic handwriting text[END_REF] 2) and 4 orientation angles (0, π/4, π/2 and 3π/4) were pre-defined. The four directions (0, π/4, π/2 and 3π/4) are widely used in the literature [START_REF] Jain | Text segmentation using Gabor filters for automatic document processing[END_REF][START_REF] Jain | Page segmentation using texture analysis[END_REF][START_REF] Zhu | Font recognition based on global texture analysis[END_REF][START_REF] Ma | Gabor filter based multi-class classifier for scanned document images[END_REF].

Designing the proper channels in order to generate filters tuned to several different frequencies and orientations has been illustrated as the crucial issue in using Gabor filters for texture characterization. Dunn et al. suggested an automatic approach for finding the optimal channels for discriminating distinct textures, but the computational complexity was very high [START_REF] Dunn | Texture segmentation using 2-D Gabor elementary functions[END_REF][START_REF] Dunn | Optimal Gabor filters for texture segmentation[END_REF]. Bianconi and Fernández evaluated the impact of the Gabor filter parameters on texture classification [START_REF] Bianconi | Evaluation of the effects of Gabor filter parameters on texture classification[END_REF]. They reported that an increase number of frequencies and orientations had an insignificant influence on texture classification performance. But, they confirmed that the best performance of texture classification is conditioned by the design of the convenient Gabor channels. Clausi and Jernigan presented a comparative study of different techniques used to extract Gabor descriptors for texture discrimination [START_REF] Clausi | Designing Gabor filters for optimal texture separability[END_REF]. The authors showed that the magnitude response outperformed the other different evaluated methods, such as using only the real component, etc. Arivazhagan et al. introduced the rotation invariant features by computing the mean and variance of the Gabor filtered image for texture classification [START_REF] Arivazhagan | Texture classification using Gabor wavelets based rotation invariant features[END_REF].

In this work, the magnitude response of the output of Gabor functions is investigated. The magnitude of the output is important if the specified Gabor filter matched the particular texture, otherwise low response to the specified Gabor filter corresponds to poor match of the dominant texture properties of the analyzed image to the set of the spatial-frequency components of the fixed Gabor filter [START_REF] Bovik | Multichannel texture analysis using localized spatial filters[END_REF]. By convoluting the analyzed whole document image at each specified channel defined by a pair of orientation and frequency, Gabor features are extracted. [START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters[END_REF] Gabor filters are applied (6 distinct spatial frequencies

f g ={2 √ 2, 4 √ 2, 8 √ 2, 16 √ 2, 32 √
2 and 64 √ 2} and 4 different orientations θ g ={0, π/4, π/2 and 3π/4}). The space of Gabor filter is set constant σ = σ x = σ y = 1. Then, two simple statistics: the mean value and standard deviation of the Gabor filtered magnitude responses corresponding to all pixels defined in the analyzed sliding window of the filtered image, are extracted. They represent the statistical distribution of Gabor magnitude responses:

• Mean: This feature characterizes the average of the Gabor filtered magnitude responses corresponding to all pixels defined in the analyzed sliding window of the filtered image. This descriptor quantifies how match the dominant texture properties of the analyzed image to the set of spatial-frequency components of the fixed Gabor filter (cf. equation ( 2)).

• Standard deviation: This descriptor determines how much the dispersion from the computed mean of the Gabor filtered magnitude response exists (cf. equation ( 3)).

The two simple first order statistics representing Gabor features are detailed in Table 1.

Features Expression

Mean F

(fg,θg) = M x=1 N y=1 I G (fg ,θg ) (x, y) M * N (1)
Standard deviation

F (2) (fg,θg) = M x=1 N y=1 [I G (fg ,θg (x, y) -F (1)
(fg,θg) ] 2 M * N (3)
In Table 1, G e(fg,θg) and G o(fg,θg) denote the spatial frequency responses of the even and odd symmetric Gabor filters.

G (fg ,θg ) (α, β) = [G e(fg ,θg ) (α, β)] 2 + [G o(fg ,θg ) (α, β)] 2 G e(fg ,θg ) = H 1(fg ,θg ) (α, β) + H 2(fg ,θg ) (α, β) 2 G o(fg ,θg ) = H 1(fg ,θg ) (α, β) + H 2(fg ,θg ) (α, β) 2j H 1(fg ,θg ) (α, β) = exp{-2πσ 2 [(α -fg cos θg) 2 + (β -fg sin θg) 2 ]} H 2(fg ,θg ) (α, β) = exp{-2πσ 2 [(α + fg cos θg) 2 + (β -fg sin θg) 2 ]} j 2 = -1

PROPOSED DOCUMENT ENHANCEMENT AND SEGMENTAION ALGORITHM

To ensure a relevant segmentation of graphical regions from textual ones on the one hand, and discrimination of text in a variety of situations of different fonts and scales on the other hand in ancient documents, a fast automatic ancient document enhancement and segmentation algorithm is proposed in this article. Figure 1 illustrates the detailed schematic block representation of the proposed algorithm.

Firstly, an ancient document image is fed as input (cf. Figure 2(a)) and is read as a gray-scale image. The different steps of the proposed algorithm is processed on gray-scale document images without introducing a binarizing task. A binarization step is avoided because it causes a loss of information specifically for information extraction and analysis. Then, to obtain enhanced backgrounds of noisy ancient documents and reduce the step complexity of a pixel-based segmentation of ancient documents, a foreground-background separation task based on SLIC segmentation [START_REF] Liu | Entropy rate superpixel segmentation[END_REF] and kmeans clustering [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] algorithms is performed.

By setting the number of superpixels ks equal to 0.1% of image pixels, the SLIC technique is carried out on a grayscale image and over-segmented image is produced. Afterwards, the background and foreground superpixels are classified based on computing the mean gray-level value of each superpixel, which is determined by averaging over all the gray-level pixels belonging to the superpixel region, and using the k-means algorithm. As a non-supervised clustering technique, the k-means algorithm is performed on the computed mean gray-level values of superpixels, without taking into account the image spatial coordinates, by setting the number of clusters kc equal to 2 to extract two clusters and subsequently to segment an image into two layers. One represents the information of the background (cf. Figure 2(c)) and the other represents the foreground (e.g. noise, text fields, drawings, etc.) (cf. Figure 2(d)). Since the k-means algorithm is carried out, the background superpixels of the original gray-level image are only processed by assigning the value of a white pixel (i.e. a 255 gray-level value) to their centers and the pixels belonging to them. However, the values of the gray-level foreground superpixels and their pixels of the original gray-level image are remained unchangeable. Thus, an enhanced and non-noisy background is achieved (cf. Figure 2(e)). Figure 2(e) illustrates an example of enhanced image by the superpixel technique with a clean background.

Afterwards, [START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters[END_REF] [START_REF] Daugman | Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by twodimensional visual cortical filters[END_REF] Gabor filters, 24 magnitude responses are generated (i.e. 24 transformed images by the selective Gabor filters are produced). Then, a quick and easy way to compute Gabor features on the whole transformed image by the selective Gabor filter, is to introduce a border replication step before the texture extraction task. By using rectangular overlapping processing windows, Gabor descriptors are only extracted from the selected foreground superpixels of the transformed image by the selective Gabor filter and the border replication step, at four different sizes of sliding windows: (16 × 16), (32 × 32), (64 × 64) and (128 × 128) to adopt a multi-scale approach. Thus, a feature vector (with dimension 48 to represent 24 Gabor filters) is produced based on the computed mean and standard deviation of the magnitude response of the transformed image by the selective Gabor filter which are extracted from one analyzed sliding window. Thus, a total of 48 features are extracted from each foreground superpixel defined in one analyzed sliding window. A 192-dimensional feature vector is subsequently formed through four different sizes of sliding windows.

Then, a foreground superpixel clustering task is performed by partitioning Gabor-based feature sets into compact and well-separated clusters in the feature space for segmenting graphical regions from textual ones on the one hand, and discriminating text in a variety of situations of different fonts and scales on the other hand. The foreground superpixel clustering task does not include spatial information and is performed by using the k-means algorithm.

Finally, a phase of labeling clusters of the gray-level foreground superpixels and the gray-level pixels belonging to each superpixel in the enhanced document image is carried out with respect to the results of the superpixel clustering phase.

Since the clustering and labeling phases of the proposed algorithm had been performed, the homogeneous regions in the ancient document are determined by labeling clusters or groups of pixels which share similar textural characteristics (cf.

EXPERIMENTS

We have experimentally evaluated the proposed ancient document enhancement and segmentation algorithm on 500 pages of ancient documents. In this section, we discuss the performance of the proposed algorithm in detail after describing our experimental corpus and its associated ground-truth.

Corpus and preparation of ground-truth

Our corpus is structured into two categories of real scanned historical documents differentiated by their content (cf.

Figure 3), reflecting the challenges of our work to determine if texture features can segment efficiently graphical regions from textual ones on the one hand, and discriminate text in a variety of situations of different fonts and scales on the other hand. Our experimental corpus includes a sufficient number of historical documents with both simple and complex layouts in each of the content categories which are ground-truthed to ensure a better understanding of the behavior of the extracted Gabor features. It is composed of:

• 250 pages containing graphics and one text font (cf. Figure 3(a))

• 250 pages containing only two fonts (cf. Figure 3(b))

The ground-truth for document images are manually outlined using rectangular regions drawn around each selected zone. The zones are ground-truthed by zoning manually each content type, i.e. each rectangular region is classified into text or graphics. Different labels for regions with different fonts are also defined for evaluating the performance of Gabor feature to separate distinct text fonts. Ground-truth is performed using the ground-truthing editor, Groundtruthing Environment for Document Images (GEDI), a public domain document image annotation tool that labels spatial boundaries of regions * [START_REF] Doermann | GEDI -A Groundtruthing Environment for Document Images[END_REF]. By specifying rectangular regions on a document image and assigning them to one of the many pre-defined content types, GEDI generates an XML schema representing the location on the page, height, width and label of each region (cf. Figure 3(c)).

Evaluation and results

The results of the proposed ancient document enhancement and segmentation algorithm using the SLIC superpixel, multi-channel Gabor filter and k-means techniques are illustrated in Figures 4(a This way of assessing the effectiveness of a segmentation method is inherently a subjective evaluation and we need to assess and evaluate robustness using an appropriate quantitative metric. Thus, in this study F-measure (F ) is computed.

To get an insight into the classification accuracy, a confusion matrix (error matrix or contingency table) [START_REF] Jensen | Introductory digital image processing[END_REF] is usually computed. From the confusion matrix, F is deduced [START_REF] Makhoul | Performance measures for information extraction[END_REF]. This accuracy metric is related to how representative the clusters are of classes and help to determine classes which are not able to segregate groups of data and give an insight into the confusion and mis-classification rates.

The computed accuracy classification values are congruent and very promising. 77% of mean F-measure is obtained by evaluating 500 pages of ancient documents. 82% and 72% of mean F-measure are noted for documents containing graphics and one text font, and only two fonts, respectively. Promising accuracy results are obtained with pages containing graphics and one text font, compared to pages containing only two fonts. So this confirmed our hypothesis that Gabor attributes have a much greater discriminant power for separating single text font and graphic regions than for separating two distinct text fonts. Nevertheless, the classification result could be improved by integrating a processing stage after the pixel labeling task, which consists of pixel grouping by taking into consideration the topographical relationships of pixels and their labels. Table 2. Evaluation of the proposed ancient document enhancement and segmentation algorithm by computing F-measure (F ). µ(F ) and σ(F ) represent the mean and standard deviation values of F , respectively when extracting and analyzing Gabor features. The " * " value is obtained by extracting and analyzing Gabor features from the selected foreground SLIC superpixels. The " * * " value is obtained by extracting and analyzing Gabor features from all foreground pixels. The experiment was run on a SGI Altix ICE 8200 cluster (1 CPU and 2 gigabytes (GB) allocated memory on a Quad-Core X5355@2.66GHz running Linux). The time required to process a page (1965 × 2750 pixels) using the proposed algorithm with extracting and analyzing Gabor features from the selected foreground SLIC superpixels, and from all foreground pixels is 04 40 and 32 36 , respectively. The amount of memory consumed when using the proposed algorithm with extracting and analyzing Gabor features from all foreground pixels is 5 times higher than with extracting and analyzing Gabor features from the selected foreground SLIC superpixels. This confirmed that the use of SLIC superpixels for analyzing Gabor features is an interesting alternative since the data cardinality is reduced and a significant gain in computation time and memory is obtained.

µ * (F ) σ * (F ) µ * * (F ) σ * * (F )
A visually comparison of the results (cf. Figure 5) when extracting and analyzing Gabor features from all foreground pixels (cf. Figure 5(b)) and from the selected foreground SLIC superpixels (cf. Figures 5(a)), indicates that better results are obtained with the case of extracting and analyzing Gabor features from the selected foreground SLIC superpixels. The outperfomance obtained when extracting and analyzing Gabor features from the selected foreground SLIC superpixels highlights our choice of using the SLIC superpixels for both document enhancement and segmentation (i.e. Gabor feature analysis). In conclusion, we observe a drop in the average of 14%(F ) when extracting and analyzing Gabor features from all foreground pixels. We have demonstrated both qualitatively and quantitatively the effectiveness of using the SLIC superpixels for both document enhancement and segmentation. Gabor features are extracted and analyzed from the selected foreground SLIC superpixels and from all foreground pixels, respectively. The pixel labeling is performed from enhanced ancient documents using the k-means algorithm by setting the number of homogeneous and similar content regions (i.e. the number of clusters) to the one defined in our ground-truth, which is equal to 2.

CONCLUSIONS AND FUTURE WORK

This article proposes an ancient document enhancement and segmentation algorithm with no hypothesis concerning the document model or the typographical parameters of the document structure. The goal of this algorithm is to find homogeneous regions within the content of digitized historical document images by extracting and analyzing texture features independently of the layout of the pages. This algorithm is based on SLIC superpixel, multi-channel Gabor filter and k-means techniques. It was evaluated on 500 pages of ancient documents with promising results. We conclude that Gabor features provide a good discrimination of the foreground layers of a document, particularly between text and graphics.

The first aspect of future work will be to use the proposed algorithm on a larger database and a wider variety of ancient documents (i.e. documents containing more than two distinct text fonts and various graphic types). This study is ongoing and will evaluate the algorithm more adequately, with more convincing experimental results. We will then study and combine statistical, geometric, model-based and other frequency texture features in order to refine the segmentation and ensure a distinction between different text fonts and various graphic types. Moreover, our results could be improved if we included fuzzy clustering methods in our clustering approach. In addition, a clustering methodology will be introduced to estimate automatically the correct number of clusters (e.g. analysis of changes in average silhouette width values computed from clusters built by using the k-means algorithm). Furthermore, by integrating a new processing stage after pixel labeling, which consists of pixel grouping that takes into consideration the topographical relationships of pixels and their labels (e.g. some operators from mathematical morphology) the pixel labeling results should be improved. This will be studied in the near future.

Figure 1 .

 1 Figure 1. Proposed ancient document enhancement and segmentation algorithm.

Figure 2 Figure 2 .

 22 Figure 2. Illustration of the intermediate results of the different steps of the proposed ancient document enhancement and segmentation algorithm: figure (a) shows an example of an ancient document image (as an input of the proposed algorithm).Figure (b) illustrates its corresponding segmented image (by using the SLIC superpixel technique) and a zoomed region of it. Figures (c) and (d) show the background and foreground SLIC superpixels, respectively. Colors assigned to the background (foreground respectively) superpixels which are illustrated in Figure (c) ((d) respectively) are randomly generated. Figure (e) depicts an enhanced image and a zoomed region of it. Figure (f) illustrates the pixel-labeled image and a zoomed region of it (as an output of the proposed algorithm (graphic regions (blue), textual regions (green)).

Figure 3 .

 3 Figure 3. Illustration of some examples of our experimental corpus and its associated ground-truth. Figure (a) illustrates an example of an ancient document containing graphics and one text font. Figure (b) shows an example of an ancient document containing only two text fonts. Figure (c) depicts an example of an associated ground-truth to an ancient printed document containing graphics and one text font.

Figure 4 .

 4 Figure 4. Examples of results of the proposed document enhancement and segmentation algorithm. The pixel labeling of the extracted Gabor features is performed from enhanced ancient documents using the k-means algorithm by setting the number of homogeneous and similar content regions (i.e. the number of clusters) to the one defined in our ground-truth, which is equal to 2. Figures (a) and (b) depict four result examples of the case when Gabor features are extracted and analyzed from the selected foreground SLIC superpixels. Figure (a) illustrates some examples of ancient document containing graphics and one text font (graphic regions (green), textual regions (blue)) while figure (b) shows some examples of ancient document containing only two fonts. Since the process is unsupervised, the colors attributed to text or graphics may differ from one document image to another.

Figure 5 .

 5 Figure 5. Examples of results of the pixel labeling of Gabor features when they are extracted and analyzed from the selected foreground SLIC superpixels and all foreground pixels. Figures (a) and (b) depict four result examples whenGabor features are extracted and analyzed from the selected foreground SLIC superpixels and from all foreground pixels, respectively. The pixel labeling is performed from enhanced ancient documents using the k-means algorithm by setting the number of homogeneous and similar content regions (i.e. the number of clusters) to the one defined in our ground-truth, which is equal to 2.

Table 1 .

 1 Gabor features.

ACKNOWLEDGMENTS

This work is supported by the ANR (French National Research Agency), under Grant ANR-10-CORD-0020, which is gratefully acknowledged. The authors would like also to thank Geneviève Cron for providing access to the Gallica digital library.


Adrien Tessier
email: adrien.tessier@inserm.fr

Julie Bertrand

Emmanuelle Marylore Chenel

Emmanuelle Comets

Comparison of nonlinear mixed effects models and non-compartmental approaches in detecting pharmacogenetic covariates

Keywords: noncompartmental analysis, nonlinear mixed effects models, penalised regression, pharmacogenetics, pharmacokinetics FWER, Family wise error rate, LOQ, Limit of quantification, NCA, Noncompartmental analysis, NLMEM, Nonlinear mixed effects model, PK, Pharmacokinetics, SNP, Single nucleotide polymorphism, α, Type I error, β, Effect size coefficient, λ, γ, ξ, Penalisation parameters

come

Comparison of nonlinear mixed effects models and non-compartmental approaches in detecting pharmacogenetic covariates

INTRODUCTION

Personalized care development should improve efficacy of drugs, and limit the risks associated with their use [START_REF] Aarons | Population pharmacokinetics: theory and practice[END_REF] , and is particularly beneficial for drugs with narrow therapeutic margin and variability in response. Pharmacogenetics [START_REF] Motulsky | Drugs and genes[END_REF][START_REF] Motulsky | Pharmacogenetics, pharmacogenomics and ecogenetics[END_REF] studies the proportion of interindividual variability in drug response that can be explained by genetic variation, investigating specifically the link between the genotype and the pharmacokinetic (PK)/ pharmacodynamic (PD) (4) phenotype. In PK/PD studies, phenotypes can be observed (e.g. residual concentrations, biomarker measurements) or estimated (e.g. total exposure or specific PK parameters).

Estimated phenotypes in PK are mainly derived through two methods. The noncompartmental analysis (NCA) [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] calculates for each subject measurements of PK exposure such as the area under the curve (AUC) or maximal and trough observed concentrations through model-free approaches. It requires rich and balanced data per individual, which makes it inappropriate for studies in advanced phases of drug development or in special populations (pediatrics…). Model-based approaches on the other hand describe the dynamic phenomena through a mathematical model and estimate primary PK parameters summarising physiological processes. These methods involve NonLinear Mixed Effects Models (NLMEM) [START_REF] Sheiner | Modelling of individual pharmacokinetics for computer-aided drug dosage[END_REF], which jointly analyse data obtained on a set of individuals to determine the typical model parameters (fixed effects), and the parameters of inter and intraindividual variability (random effects), as well as the residual error. NLMEM are better suited to sparse and unbalanced data, and clinical studies can be combined to increase the power to detect genetic effect [START_REF] Sheiner | Pharmacokinetic/pharmacodynamic modeling in drug development[END_REF].

Unlike monogenic diseases, caused by mutation of a single gene, the variability in PK/ PD is usually the product of a set of markers with low to intermediate effect sizes. The screening of a large number of genetic markers, such as Single Nucleotide Polymorphisms (SNP), is possible thanks to the development of genotyping methods and microarrays [START_REF] Shalon | A DNA microarray system for analyzing complex DNA samples using two-color fluorescent probe hybridization[END_REF]. This has two major consequences: the number of genetic variants studied tends to be larger than the number of subjects and some of these variants are correlated due to linkage disequilibrium [START_REF] Daly | High-resolution haplotype structure in the human genome[END_REF].

An exhaustive literature review about clinical pharmacogenetic studies was performed to gain an overview of the methods currently used in this type of study (see the supplementary materials for the literature review methods we used). Information about the phenotypes (obtained by NCA or modelling), the design of studies, the genetic data and the statistical analysis have been extracted from each publication, and summarised through descriptive statistics. During the period 2010-2012, on 85 pharmacogenetic studies using PK parameters as phenotypes, 69% used NCA and 31% used modelling-based phenotypes for association analysis. About two thirds of the studies included less than 50 subjects and only 15% more than 100 subjects. A limited number of genetic covariates were studied (3 genes (minimum = 1; maximum = 45) and 8 SNPs (minimum = 1; maximum = 198) in average). Genes were finely targeted because known to be involved in the PK of studied molecules. Finally, association between NCA-based phenotype and polymorphism were mainly investigated using univariate methods (univariate anova, t-test…) (57%), multivariate methods (multivariate linear regression) (14%) and stepwise or descriptive methods. Model-based phenotypes were explored using mostly stepwise regression (78%), univariate methods applied on individual parameter estimates (7%) or descriptive methods.

So although health authorities strongly recommend studying the pharmacogenetics of new chemical entities in development [START_REF]Guideline on the use of pharmacogenetic methodologies in the pharmacokinetic evaluation of medicinal products[END_REF][START_REF] Fda | Guidance for Industry and FDA Staff: Pharmacogenetic Tests and Genetic Tests for Heritable Markers[END_REF], there is no consensus on analysis methods to explore a large number of polymorphisms in association with PK phenotypes. Specific approaches and statistical tools are required, which must take into account the small amount of PK information provided by each individual in relation to the number of genetic covariates, and be able to detect a signal in a large number of possible relationships.

Penalised regression methods analyse all markers simultaneously and use a penalisation function which shrinks most effect coefficients to select a parsimonious set of markers of the phenotype variability. These methods have been especially used in Genome Wide Association Studies (GWAS) to binary outcomes [START_REF] Omoyinmi | Association of the IL-10 gene family locus on chromosome 1 with juvenile idiopathic arthritis (JIA)[END_REF] or quantitative traits [START_REF] Yao | Genome-wide association study of lung function phenotypes in a founder population[END_REF]. Here we evaluate three such methods: ridge regression ( 14) adapted to include a test of significance for its derived estimates [START_REF] Cule | Significance testing in ridge regression for genetic data[END_REF], the Lasso ( 16) and the HyperLasso [START_REF] Hoggart | Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies[END_REF], a generalisation of the Lasso.

In the present study, we propose to compare those three penalised regression methods with a stepwise approach through a simulation study, to assess their ability to detect the influence of genetic variables on the PK. We apply them to four possible PK phenotypes (two observed concentrations, AUC estimated by NCA, model parameters estimated by NLMEM). This work is based on collaboration with the pharmaceutical industry (IRIS, Institut de Recherches Internationales Servier) and thus focuses specifically on issues of clinical trials for drug development. This also provided a real case study to design the simulations which had interesting features, including a nonlinear absorption, resulting from the drug's physicochemical properties, and a large number of SNPs collected thanks to a specific microarray. This could enable a meaningful comparison of phenotypes in a challenging context.

MATERIALS AND METHODS

Overview

In this work, we use a model previously developed to fit the data from a real example for a molecule in early clinical development at IRIS (Motivating data example section). For this work we have simulated PK profiles using the model and parameters based on the real data to derive phenotypes (PK phenotypes section), on which we have applied different association methods (Statistical methods for genetic association section). We have simulated scenarios (Simulation study section) for different experimental protocols varying especially the number of subjects and evaluate the methods in term of detection probability (Evaluation section).

Motivating data example

The simulation study was inspired by a real dataset which we cannot use for confidentiality issues, but which presented characteristic features encountered during the early phases of drug development, including frequent PK sampling, extensive genetic data, a complex PK profile and a wide range of doses. We start by describing the clinical study and the experimental protocols, as well as the model that will be used in the simulation study.

Drug S developed by IRIS was investigated in three phase I clinical studies performed in a total of 78 adult healthy volunteers. All subjects were genotyped at baseline using a DNA microarray developed by the laboratory. This chip has been developed specifically for PK studies and consists only of SNPs known for being involved in the PK of drugs, i.e. markers of phase I and II metabolic enzymes, of SLC or ABC family transports, as well as nuclear receptor genes. The studies included 8 dose groups [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF][START_REF]Guideline on the use of pharmacogenetic methodologies in the pharmacokinetic evaluation of medicinal products[END_REF][START_REF] Duffull | Evaluation of the pre-posterior distribution of optimized sampling times for the design of pharmacokinetic studies[END_REF]50,100,200,400 or 800 units) with different allocations (respectively 6, 6, 24, 12, 12, 6, 6 and 6 subjects by dose) and extensive PK sampling.

A two-compartment model with a double absorption function (Fig. 1, Table I) was used to describe the PK of drug S. The nonlinearity with dose was due to the low solubility of drug S, and was modelled through a dependency of the absorption parameters F and FRAC on dose.

The relationship of F with dose was expressed through an Imax model [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] parameterised in Imax F and D50 F while the relationship of FRAC with dose was expressed with an Emax model [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] parametrised in Emax FRAC and D50 FRAC . Tk0 is the zero order absorption constant rate and

Pharmacokinetic phenotypes

To capture the PK signal, we considered 3 types of measures in this study.

Observed concentrations

Observed concentrations are easy to obtain and require no additional analysis step. Here we considered the last concentration at 192h after a single dose administration (C192h), and the concentration at 24h (C24h), which corresponds to a trough concentration when repeated doses are given during routine treatment.

Noncompartmental approach

AUC was calculated using the linear trapezoidal rule and extrapolation to infinity was achieved assuming an exponential decay [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF][START_REF] Jaki | Estimation of pharmacokinetic parameters with the R package PK[END_REF].

Normalisation of observed and NCA phenotypes

To take into account the nonlinearity in dose absorption, observed concentrations C24h and C192h, and AUC obtained by NCA were normalised by dose. An Emax model [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] where 𝜃 𝑖 is the vector of individual parameters assumed to follow a log-normal distribution:

𝜃 𝑖 = 𝜇𝑒 𝜂 𝑖
with 𝜇 the population average parameter and 𝜂 𝑖 the difference between the population average and the individual 𝑖 parameter. We assumed that PK parameters follow a log-normal distribution to ensure that they are strictly positives [START_REF] Duffull | Evaluation of the pre-posterior distribution of optimized sampling times for the design of pharmacokinetic studies[END_REF].

𝜀 𝑖𝑗 is the residual effect that quantifies the deviation between the model prediction and the measured value for subject 𝑖 at time𝑗. The residual variability was described using a proportional model:

𝑔(𝜃 𝑖 ; 𝑡 𝑖𝑗 ) = 𝜎 𝑠𝑙𝑜𝑝𝑒 (𝑓(𝜃 𝑖 ; 𝑡 𝑖𝑗 ))
where 𝜎 𝑠𝑙𝑜𝑝𝑒 represents the proportional component. We typically assume 𝜂 𝑖 ~𝑁(0, 𝜔 2 ) and 𝜀 𝑖𝑗 ~𝑁(0, 𝜎 2 ).

Statistical methods for genetic association

We assume that a linear model links the phenotype to the genetic variants as in:

𝑃ℎ𝑒𝑛𝑜𝑡𝑦𝑝𝑒 𝑖 = 𝛽 0 + ∑ 𝛽 𝑘 . 𝑆𝑁𝑃 𝑖𝑘 + 𝜀 𝑖 ; 𝑆𝑁𝑃 𝑖𝑘 = {0, 1, 2}
where 𝑃ℎ𝑒𝑛𝑜𝑡𝑦𝑝𝑒 𝑖 is a vector for individual phenotypes, 𝛽 0 an intercept, 𝑆𝑁𝑃 𝑖𝑘 a vector for the genetic variants and 𝛽 𝑘 a vector for the individual genetic effect size associated to the genetic variant and 𝜀 𝑖 a residual error following a Gaussian distribution. In this model, the genetic variant takes values 0, 1 or 2, reflecting the number of mutated alleles. All phenotypes are log-transformed to ensure that they follow a normal distribution. To account for type I error inflation due to the multiplicity of tests, all methods used a Sidak correction on the Family Wise Error Rate (FWER) to compute a type I error per SNP α:

𝛼 = 1 -(1 -𝐹𝑊𝐸𝑅) 1 𝑁 𝑡 ×𝑃 𝑡
where 𝑁 𝑡 and 𝑃 𝑡 are the numbers of SNPs and PK parameters considered simultaneously.

Ridge regression

Ridge regression imposes a penalty on the size of the 𝛽 𝑘 to reduce the prediction error [START_REF] Cule | Significance testing in ridge regression for genetic data[END_REF] without preventing the integration of the model variables. We used the approach proposed by Cule et al. to set semi-automatically the penalty so that the trace of the projection matrix, which relates the predictions to the observations, is equal to the number of components in a principal component analysis (PCA) of the data. From a Bayesian perspective, this correspond to applying a Gaussian prior of identical variance on the eigenvalues issued from the PCA [START_REF] Cule | Ridge regression in prediction problems: automatic choice of the ridge parameter[END_REF]. The ridge regression do not shrunk coefficients estimates to 0. Therefore we used a Wald test on these coefficients and their standard error (SE), as proposed by Cule et al. [START_REF] Cule | Significance testing in ridge regression for genetic data[END_REF] to perform the variable selection with the test statistic 𝑇 0 :

𝑇 0 = 𝛽 𝑘 𝑆𝐸(𝛽 𝑘 )
where 𝑆𝐸(𝛽 𝑘 ) is the standard error of the regression coefficient 𝛽 𝑘 . Under the null hypothesis, 𝑇 0 follows a Student t distribution with a significant threshold equal to the type I error per SNP α.

Lasso

Lasso ( 16) also uses a penalty function, which from a Bayesian point of view corresponds to using a double-exponential (DE) probability density as a prior on 𝛽 𝑘 . The Lasso sets some coefficients to 0 for sufficiently large values of the tuning parameter; this allows variable selection and ensures a parsimonious model. The regularisation parameter 𝜉 is calculated to achieve a target FWER, using the following expression [START_REF] Hoggart | Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies[END_REF]:

𝜉 = 𝛷 -1 (1 - 𝛼 2 ) √ 𝑁 𝜎 2
where α is the type I error per SNP, 𝛷 -1 is the inverse normal distribution function, N the number of subjects and σ the standard error of the phenotype considered.

HyperLasso

HyperLasso ( 17) is derived from the Lasso. Here the penalisation corresponds to using a normal-exponential gamma (NEG) distribution as a prior on 𝛽 𝑘 and depends on two parameters: a shape parameter 𝜆 and a scale parameter 𝛾. The sharp peak at zero and the flatter tail of the NEG distribution favour sparse solutions but the estimates of larger effects are shrunken less severely than the Lasso. The smaller the shape parameter the heavier the tails of the distribution and the more peaked at zero, which can result in fewer correlated SNPs being selected. The shape parameter 𝜆 was set to 1 in our study, which gives realistic effect size distributions [START_REF] Vignal | Using penalised logistic regression to fine map HLA variants for rheumatoid arthritis[END_REF]. As for Lasso, the scale parameter 𝛾 is calculated, again depending on a type I error per SNP α, using the following expression [START_REF] Hoggart | Simultaneous analysis of all SNPs in genome-wide and re-sequencing association studies[END_REF]:

𝑠𝑖𝑔𝑛(𝛽)(2𝜆 + 1) 𝛾 𝐷 -(2𝜆+2)( |𝛽| 𝛾 ) 𝐷 -(2𝜆+1)( |𝛽| 𝛾 ) = 𝛷 -1 (1 - 𝛼 2 ) √ 𝑁 𝜎 2
where α is the type I error per SNP and D the parabolic cylinder function [START_REF] Gradshteyn | Tables of Integrals, Series and Products: Corrected and Enlarged Edition[END_REF].

Stepwise procedure

We use the algorithm in Figure 2 inspired by Lehr et al. [START_REF] Lehr | Integration of high-throughput genotyping data into pharmacometric analyses using nonlinear mixed effects modeling[END_REF] and based on univariate regression: (i) PK phenotypes are regressed on each SNP and a Wald test is applied with a significance threshold equal to α. To account for linkage disequilibrium, among selected SNPs showing strong correlation (r² > 0.8), only the most significant is kept. Then (ii) the most significant among selected SNPs is included in the linear regression of the PK phenotype on the SNPs. These iterative steps are performed until no more SNP enter the linear model.

Simulation study

Genotypes

SNPs were simulated using the Hapgen2 software [START_REF] Su | HAPGEN2: simulation of multiple disease SNPs[END_REF] based on the DNA microarray used in clinical studies. To simulate genetic variants for the 200 data sets while retaining the correlations between variants found in the human genome, we used a reference panel of Hapmap genotypes data set (Hapmap 3 release 2) for a Caucasian population [START_REF]The International HapMap Project[END_REF]. Hapgen2

simulates genotype with the same LD patterns as the reference data [START_REF] Su | HAPGEN2: simulation of multiple disease SNPs[END_REF].

Summaries of simulated genotypes, i.e. minor allele frequencies, Hardy-Weinberg equilibrium and LD plots, can be found in supplementary file (Simulated polymorphisms information).

Pharmacokinetic profiles

Using the model and parameters estimates from the motivating data example (Table I), we simulated concentrations after a single dose administration according to an extensive sampling schedule with 16 samples per subject at 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 16, 24, 48, 72, 96, 120 and 192 hours after taking the tablet. We did not include a limit of quantification (LOQ) and the data were simulated without censoring. NCA method was applied on these simulated profiles and for the model-based approach, population and individual parameters were re-estimated using the Monolix software [START_REF] Lavielle | The MONOLIX software[END_REF] and SAEM estimation algorithm [START_REF] Kuhn | Coupling a stochastic approximation version of EM with an MCMC procedure[END_REF].

Individual CL, Q and V2 were used as PK phenotype in a post-model covariate analysis step.

In each simulation scenario, we performed two simulations: in the first simulation (null hypothesis H 0 ) we assumed there was no effect of genetics on the PK parameters; in the second simulation (alternative hypothesis H 1 ), 6 SNPs were drawn randomly and we assumed they had an impact on the log-transformed clearance according to the following model:

𝑙𝑜𝑔(𝐶𝐿𝑖) = 𝑙𝑜𝑔(𝜇 𝐶𝐿 ) + ∑ 𝛽 𝑘 × 𝑆𝑁𝑃 𝑖𝑘 6 𝑘=1 + 𝜂 𝑖 𝐶𝐿
where 𝐶𝐿𝑖 is the individual clearance, 𝜇 𝐶𝐿 the population clearance, 𝛽 𝑘 the effect size associated to the genotype 𝑆𝑁𝑃 𝑖𝑘 and 𝜂 𝑖 𝐶𝐿 the interindividual variability in clearance.

For each SNP, the associated effect size 𝛽 𝑘 is computed as a function of the coefficient of genetic component (𝑅 𝐺𝐶 𝑘 ) and the minor allele fraction (𝑝 𝑘 ), according to the following equation:

𝛽 𝑘 = √ 𝑅 𝐺𝐶 𝑘 × 𝜔 2 𝐶𝐿 2𝑝 𝑘 (1 -𝑝 𝑘 ) -𝑅 𝐺𝐶 𝑘 × 2𝑝 𝑘 (1 -𝑝 𝑘 )
where 𝑅 𝐺𝐶 𝑘 is the part of the interindividual variability in CL explained by the SNP (expressed in %) and 𝜔 2 𝐶𝐿 is the variance of random effects on CL due to non-genetic sources. To simulate realistic genetic effects, the 6 SNPs altogether explained a total 𝑅 𝐺𝐶 of 30% with unbalanced effect sizes (𝑅 𝐺𝐶 𝑘 respectively equal to 1, 2, 3, 5, 7 and 12%). These effect sizes were chosen to be consistent with genetic effects observed in clinical studies. For example warfarin doses has been associated with three genetic variants in the cytochrome P450 warfarin-metabolizing genes CYP4F2 and CYP2C9 and in the warfarin drug target VKORC1, explaining respectively 1.5, 12 and 30% of variability [START_REF] Takeuchi | A genomewide association study confirms VKORC1, CYP2C9, and CYP4F2 as principal genetic determinants of warfarin dose[END_REF].

Simulation scenarios

We simulated different scenarios. In the first scenario S real , the design was chosen to be close to the drug S Phase I clinical trials protocol: 78 subjects (N) receiving 8 different doses [START_REF] Gabrielson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF][START_REF]Guideline on the use of pharmacogenetic methodologies in the pharmacokinetic evaluation of medicinal products[END_REF][START_REF] Duffull | Evaluation of the pre-posterior distribution of optimized sampling times for the design of pharmacokinetic studies[END_REF]50,100,200, 400 or 800 units) with a similar allocation ratio and 16 sampled times (n) at 0.5, 1, 1.5, 2, 3, 4, 6, 8, 12, 16, 24, 48, 72, 96, 120 and 192h. A second scenario S large using the same design but with more subjects (N = 384, n = 16) was considered as well to investigate a large sample situation. This scenario represents an ideal case from a design perspective, with a large number of subjects and a rich protocol for each subject. Figure 3 To evaluate the effect of the PK model structure including the nonlinearity of absorption, 3 additional scenarios were simulated using the same design than S large but different structural models: a nonlinear absorption without interindividual variability on F (S large, no IIV_F ), a linear PK (S linearPK : FRAC=0.5, F=0.8, (F)=32.9%) and a linear PK without variability on F (S linearPK, no IIV_F : FRAC=0.5, F=0.8, (F)=0).

Evaluation

Each method: ridge regression, Lasso, HyperLasso and the stepwise procedure was applied to each PK phenotype: observed C24h and C192h, AUC estimated by NCA and model parameters CL, V and Q estimated by NLMEM. To enable proper power comparison the target FWER was set to 20% (with a prediction interval for 200 data sets equal to [14.5-25.5]). Under H 0 , an empirical FWER was estimated as the percentage of data sets where at least one SNP was found significant. If necessary, the type I error per SNP α were corrected empirically until the FWER estimate was not significantly different from 20%. Under H 1 , we recorded for each method and phenotype the number of true positives (TP, corresponding to the selection of a SNP which was indeed associated to CL in the simulation, its maximum over the 200 simulations being 1200) and false positives (FP, corresponding to a SNP selected in the model but not present in the simulation). A 95% confidence interval was also estimated assuming the number of TP or FP follows a Poisson distribution. The true positive rate (TPR) and the false positive rate (FPR) were calculated as follow:

𝑇𝑃𝑅 = 𝑇𝑃 𝑇𝑃 + 𝐹𝑁 𝐹𝑃𝑅 = 𝐹𝑃 𝐹𝑃 + 𝑇𝑁
where FN is the count of false negatives and TN the count of true negatives. The TPR was the main outcome on which we statistically compared the different methods and phenotypes.

First, we compared the TPR across the different phenotypes for each method, using Cochran's Q test [START_REF] Cochran | The Comparison of Percentages in Matched Samples[END_REF] in the following way: for each dataset and each phenotype, we defined a variable X with value 1 for the phenotype(s) with the maximal TPR, and 0 for the other ones; this was done separately for each method. A global test for each method was performed first, with a significance threshold set to 5%; if the phenotypes were found to be significantly different, pairwise comparisons were then performed with a Wilcoxon test, using a Sidak correction for the number of tests (a corrected threshold of 0.009 for 6 pairwise comparisons).

Once the phenotype yielding the highest TPR across the different methods has been selected, the same approach was applied to compare methods on this phenotype. In each simulation, we also recorded the number of datasets for which all methods detected the same number of SNPs (X=1 for all methods), and the number of datasets for which all methods failed to detect at least one SNP (X=0 for all methods).

RESULTS

Scenario 1: S real

Control of FWER under H 0

Table II shows the estimates of the empirical FWER under H 0 . All methods tended to be too conservative since the FWER was lower than expected for all methods and phenotypes, except on the last concentration parameter C192h. After an empirical correction of thresholds or penalisation parameters depending on the method, FWER was properly controlled around 20%, as shown in Table II. This correction was applied in the corresponding simulation under H 1 .

To investigate why the FWER was not controlled under H 0 despite the calibration step within each method, we simulated a scenario S independent similar in design to S real but with independent SNPs. The FWER estimates in S independent were non-significantly different from 20% (results in supplementary materials). This suggests that correlations between SNPs lead to a decrease in FWER.

Performance under H 1

In each scenario, the total number of SNP associated with the PK was 1200, but only a small number was effectively detected, as shown in Table III. This was even more apparent for associations based on NCA parameter or observed concentrations. The TPR was significantly different between the different phenotypes (p<0.001 for each of the 4 methods, according to Cochran's Q test) and was highest for CL in each 2 by 2 phenotype tests (p<0.002 for all pairwise comparisons). The FPR on the other hand was similar for all phenotypes (figure 4a).

The methods were then compared for CL, and the TPR was found to be significantly different between the four methods (p<0.001). Ridge regression yielded a higher TPR more often than other methods (p<0.003 for all pairwise comparison), while Lasso, HyperLasso and the stepwise procedure were comparable.

The probability to detect at least one genetic variant on CL estimates was low, around 40% for all methods (Fig. 5a). This probability decreased quickly when trying to detect more variants and reached 0 for 3 variants or more.

Scenario 2: S large

Control of FWER under H 0

As expected due to the larger number of subjects, the estimated FWER increased with comparison to the scenario S real , but remained below the target of 20% for some methods.

Empirically corrected thresholds and penalty terms were determined for each combination of phenotype and method to obtain FWER estimates of 20% and were then used for the simulations under H 1 (results in supplementary materials).

Performance under H 1

The number of TP increased clearly compared to the previous scenario (Table IV).

Concerning the comparison between phenotypes we found similar results as for S real : all methods were more powerful when applied to CL compared to the other phenotypes (p<0.001 for all pairwise comparison). The number of FP also increased and was quite high, but increased proportionally less than the number of TP for the model-based phenotypes, C192h and AUC. For C24h on the other hand, the number of FP exceeded the number of TP.

Thus for a similar FPR between all phenotypes, the TPR was higher for CL (Fig. 4b). The TPR for ridge regression was higher than Lasso and HyperLasso (p<0.001 with methods applied to CL). The TPR for the stepwise procedure was intermediate, the 2 by 2 tests were not significant for the comparison of the stepwise procedure versus ridge regression or Lasso.

As expected, with increasing the number of subjects the power of all methods increased to reach almost 100% to detect at least one genetic variant (Fig. 5b). Then, the power decreased when trying to detect more variants. Departure in methods was observed on the power to detect at least 3 and more variants, with the ridge regression and stepwise procedure showing higher power.

Influence of the structural PK model under H 1

In both scenarios, regardless of the association method, the power to detect a gene effect was higher using PK parameter obtained by NLMEM than AUC estimated by NCA or observations. We investigated whether this was due to the specific features in the PK model, i.e. the non-linearity in the absorption model and/or the variability in the bioavailability.

When we assumed no variability on the bioavailability parameter F (scenario S large, noIIV_F ), the difference in the number of TP between CL and AUC was smaller (Fig. 4c), but the TPR for CL remained higher compared to AUC. Assuming a linear absorption, while retaining variability on F, also reduced this difference (scenario S linearPK , Fig. 4d) but only when we used a linear absorption model without variability on F did the benefit of CL over AUC disappear (scenario S linearPK, noIIV_F , results in supplementary materials). Changes in the PK structural model did not affect the number of TP for C192h, while removing the variability on bioavailability increased the number of TP for C24h, although it remained very low.

DISCUSSION

Many analysis methods have been proposed in the pharmacogenetic literature, depending on the phenotypes studied and the association test. NCA is mainly used to represent the PK exposure, univariate association methods are still widely applied and sample sizes are limited. This work aimed to evaluate 16 combinations of four phenotypes and four methods and the design of the different simulation scenarios was based on actual clinical studies. This realistic setting enabled a meaningful comparison, providing at the same time a challenging context of nonlinear PK and a custom set of polymorphisms. This work takes place in context of exploratory analyses, and we therefore chose a high FWER for variable selection (20%). The four association methods use estimated phenotypes after an initial estimation step without covariates included in the model. Work on alternative approaches which simultaneously estimate the PK model parameters and the genetic size effects are ongoing [START_REF] Bertrand | Penalized regression implementation within the SAEM algorithm to advance high-throughput personalized drug therapy. 22th PAGE Meet Glasg Scotl[END_REF]; they require iterative selection and estimation which increases the computational burden. We decided to study the most common association methods based on a maximum likelihood approach, the ridge regression and Lasso, together with a specific extension for genetic covariates (HyperLasso). Other penalized regression methods have been proposed such as the elastic net method [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], which has shown intermediate performances between ridge regression and Lasso. Other approaches to investigate include modifications of the penalised methods we used, such as the significance test very recently developed for the Lasso [START_REF] Lockhart | A significance test for the lasso[END_REF]. In the present work, we did not consider gene-gene and geneenvironment interactions. Model-based approaches have been proposed in such contexts, and evaluated on real pharmacogenetic data sets [START_REF] Knights | Vertical Integration of Pharmacogenetics in Population PK/PD Modeling: A Novel Information Theoretic Method[END_REF]. These methods or clustering based algorithms should be compared to penalized regression methods in simulations close to those presented in this work. Frequentist approaches are most often used in population PK/PD analysis. Bayesian methods [START_REF] O'hara | A review of Bayesian variable selection methods: what, how and which[END_REF] may also be worth considering for variable selection but there use remains limited in population PK/PD studies. The estimate of entire distributions of parameters adds an additional numerical complexity requiring further development beyond the scope of the present work.

Several methodological works associating pharmacogenetics and NLMEM have been published [START_REF] Bertrand | Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters[END_REF][START_REF] Bertrand | Pharmacogenetics and population pharmacokinetics: impact of the design on three tests using the SAEM algorithm[END_REF][START_REF] Bertrand | Some alternatives to asymptotic tests for the analysis of pharmacogenetic data using nonlinear mixed effects models[END_REF][START_REF] Bertrand | Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models[END_REF]. Lehr et al. have suggested an adaptation of the classical stepwise covariate selection on PK phenotype in NLMEM [START_REF] Lehr | Integration of high-throughput genotyping data into pharmacometric analyses using nonlinear mixed effects modeling[END_REF] and a method inspired from Lasso has already been used for the selection of non-genetic covariates in NLMEM [START_REF] Ribbing | The lasso--a novel method for predictive covariate model building in nonlinear mixed effects models[END_REF]. But this is the first work comparing model-based approach with NCA in this area. In our study, all methods showed a higher number of TP when used on individual clearances CL from NLMEM, compared to the other phenotypes (AUC, C24h, and C192h). Furthermore, relatively to the number of TP, the number of FP was lower for model-based phenotypes (CL, Q and V2) than other phenotypes, also improving the TPR. This finding indicates that using a modelling approach enables better power. Indeed, the modelling approach allows separating the different phases of the PK process (absorption, distribution, metabolism, excretion), improving the interpretation of the genetic effects by the comprehension of mechanisms behind the association between a genetic and a primary PK parameter. The benefits provided by this approach compared to the NCA in terms of power have also been shown in other areas [START_REF] Dubois | Bioequivalence tests based on individual estimates using non-compartmental or model-based analyses: evaluation of estimates of sample means and type I error for different designs[END_REF][START_REF] Panhard | Evaluation by simulation of tests based on non-linear mixedeffects models in pharmacokinetic interaction and bioequivalence cross-over trials[END_REF], especially when the number of samples per subject is limited, but remained to be demonstrated in the field of pharmacogenetics. In the simulation we did not include a LOQ and the data were simulated without censoring. In practice however late measurement of last concentrations could entail a significant number of data below LOQ which would decrease the ability to detect genetic differences. We would expect the NCA approach to be also impacted, as data below the LOQ are usually omitted in NCA, resulting in bias in parameters estimated through such approach [START_REF] Humbert | Evaluation of pharmacokinetic studies: is it useful to take into account concentrations below the limit of quantification?[END_REF]. In NLMEM on the other hand different methods based on likelihood have been developed to impute values below the LOQ in NLMEM [START_REF] Ahn | Likelihood based approaches to handling data below the quantification limit using NONMEM VI[END_REF], resulting in unbiased estimates. With such methods, presence of data below the LOQ should not modify the probability to detect genetic effects on phenotypes estimated through NLMEM.

AUC is highly correlated with CL (𝐴𝑈𝐶 = 𝐷𝑜𝑠𝑒 𝐶𝐿 ) when the number of samples per subject is large, as in design used in simulations. Despite this the power to detect a gene effect was higher for the model-based approach than for the phenotypes estimated by NCA or observed due to the specific features in the PK model, the nonlinearity in the absorption model and the variability in the bioavailability. Indeed, with a linear absorption and no variability on the bioavailability parameter F (scenario S linearPK, noIIV_F ), the TPR is similar between CL and AUC (results in supplementary materials). Still in the linear case, an interindividual variability on F reduces the TPR of AUC (scenario S linearPK ). In the nonlinear case, although the phenotypes observed or estimated by NCA are normalised, because of estimation errors in the Emax model used for the normalisation, their respective TPR are lower than CL, even with no variability on F (scenario S large, noIIV_F ). In such a rich design, where the subjects are extensively sampled, AUC is appropriate when the PK is simple (linear, without variability on the absorption process).

In our simulations we introduced effects from six SNPs on CL parameter in the PK model to simulate concentrations profiles under the alternative hypothesis. These settings may favour the model-based phenotypes and also by extension AUC which is directly correlated with CL.

With a nonlinear PK model model-based phenotypes proved much more powerful compared to the phenotype estimated by NCA, even correcting the AUC by dose to take into account the nonlinearity. On the other hand, when the PK model was linear, we found similar results in terms of TPR and FPR for both phenotypes, but only when there was no interindividual variability on bioavailability. Concerning the observed phenotypes we expected the late concentration C192h to be a good reflection of clearance, and to give similar performance than CL after correction for nonlinearity. However, this was not the case in our results, even with a linear PK. The reason for this is probably that the influence of the other parameters dilutes the impact of the genotypes on this observed phenotype, while AUC only depends on CL. But the results in term of TP obtain on C192h were better than on C24h. This concentration is not informative for elimination clearance because occurring for many subjects during the rebound simulated using our PK model, so that the poor performance for this phenotype could be expected.

In the model-based approach, using NLMEM individual parameter estimates are derived after the estimation of the population parameters and their precision depends on the amount of individual information available in the data [START_REF] Savic | Importance of shrinkage in empirical bayes estimates for diagnostics: problems and solutions[END_REF]. In the scenarios we simulated, the number of sampling points per subject was large, so that all model parameters, included the phenotype, were well estimated. The number of TP obtained using the simulated clearances (without the estimation step) in both scenarios was only slightly higher than when using the estimated clearance.

The estimates of the probability to detect causal variants were similar between methods in scenario S real , while a departure was noticed in scenario S large to detect at least 3 variants.

The ridge regression method exhibited the highest count of TP but to the cost of a higher number of FP. In scenario S real the estimates were rather low and fell to nearly 0% for the probability to detect at least 3 of the 6 SNPs, indicating a low ability to detect multiple effects in a realistic design. As expected, increasing the number of subjects (scenario S large ) strongly increased the probability of the different methods to detect more variants. But detecting all the 6 simulated variants remained very rare, due to the small percentage of variability explained by some variants. Our simulations illustrate the importance of the number of subjects in pharmacogenetic studies : infrequent mutations are unobserved in a small population and association methods are sensitive to the frequency of the variant allele [START_REF] Bertrand | Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models[END_REF].

All methods were fast to run, requiring several minutes to several hours depending on the number of subjects in the scenario. Despite its iterative algorithm, stepwise procedure using univariate regression was the simplest and fastest method in all scenarios. However, the use of a full stepwise procedure using model estimation like proposed by Lehr et al. [START_REF] Lehr | Integration of high-throughput genotyping data into pharmacometric analyses using nonlinear mixed effects modeling[END_REF] causes a sharp increase in run time [START_REF] Bertrand | Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models[END_REF]. Ridge regression and Lasso were intermediate methods in term of run time and HyperLasso was the slowest method.

There have been few papers comparing association methods in pharmacogenetics applied to PK/PD. The present work complements a previous study by Bertrand and Balding (39), who compared four association methods (ridge regression, Lasso, HyperLasso and a stepwise procedure) on only one kind of PK phenotype, CL estimated by NLMEM. Their setup is close to our scenario S large in terms of the number of subjects (N=300), but each subject was sampled only 6 times. We compared the number of TP and FP between our studies. Our results for model-based phenotypes are partly different from those presented in this study.

Bertrand and Balding found both fewer TP and less FP than in the S large scenario. These differences result from our respective calculations of TP and FP. In Bertrand and Balding [START_REF] Bertrand | Multiple single nucleotide polymorphism analysis using penalized regression in nonlinear mixed-effect pharmacokinetic models[END_REF] causal variants were removed from the analysis data sets and TP were defined as SNPs correlated with the causal variants with an r²>0.05. In our study on the other hand, the causal variants were present in the analysis data sets. Using the calculations of Bertrand and Balding, we obtained similar numbers of FP but the numbers of TP remained higher.

Bertrand and Balding also explored 1227 SNPs, i.e. 7 times more SNPs in a similar number of subjects. In the ridge regression, the threshold for the Wald test proposed by Cule et al. [START_REF] Cule | Significance testing in ridge regression for genetic data[END_REF] was therefore much more stringent which could explain that ridge regression detected fewer TP than the other shrinkage based approaches in their study.

In conclusion this work shows the critical importance of using modeling approaches for pharmacogenetic studies. They allow detecting associations between genetic variants and PK more efficiently, in particular in the presence of complex PK involving non-linearity, where the AUC even corrected for the dose effect was much less sensitive to the genetic effect. In addition the use of models allows for the analysis of intrinsic parameters with physiological meaning as an elimination or an absorption rate. Our results also reinforce the importance of the number of subjects in pharmacogenetic studies, and suggest that it may not be reasonable to expect to detect even strong genetic effects and/or genetic effect due to rare alleles with small sample sizes. In addition, this work has highlighted statistical difficulties; FWER was not properly controlled and lower than expected. An empirical correction had to be performed to target 20% for the FWER. This decrease was due to the correlation between SNPs, as we showed in an additional simulation with uncorrelated genotypes. Consequently, to enable the comparison across the different methods, the FWER was set empirically for the 16 combinations of methods and phenotypes . A final message from this work is that in our simulation settings no association method showed a much higher power than the others. Stepwise procedure after empirical correction b 18.5 23.5 19 18 a. Set of empirical family wise error rates (FWER) obtain without correction. b. Set of empirical FWER obtain after correction of thresholds or penalisation parameters. The 95% prediction interval around 20 for 200 simulated data sets is [14.5-25.5].

Tlag 1

 1 the corresponding lag time, Ka the first order absorption constant rate and Tlag 2 the corresponding lag time, V1 the central compartment volume, V2 the peripheral compartment volume, Q the intercompartmental clearance and CL the elimination clearance. Interindividual variability was described by an exponential model on F, Tk0, Tlag2, V2, Q and CL.

  represents the PK profiles from a simulated data set with S large design as function of the genotypes of the SNP explaining 12% of the clearance interindividual variability (under the alternative hypothesis H 1 ). The profiles show an important interindividual variability. As expected last concentrations are lower in heterozygotes and rare homozygotes under H 1 , due to the increase of clearances.

  The probability to detect a given number x (x=1, …, 6) of the 6 SNPs which were associated to CL was computed as the percentage of data sets simulated under H 1 where x or more SNPs are selected. Of note, for the model-based analysis, associations were explored on parameters CL, Q and V2. TP were causal variants associated to CL and FP were non-causal variants associated to CL and any variants associated to Q and V2. The FPR for the modelbased phenotype designed by CL were then computed by taking into account the FP on any of the three model parameters.

Table I .

 I Population values (µ) and interindividual variability (ω) for drug S model (units not disclosed for confidentiality reasons) 𝐹𝑅𝐴𝐶 = 𝐸𝑚𝑎𝑥 𝐹𝑅𝐴𝐶 𝑑𝑜𝑠𝑒 𝐷50 𝐹𝑅𝐴𝐶 + 𝑑𝑜𝑠𝑒 , where dose is the administrated amount.

	Parameters		µ	ω (%)
	F a	D50 F	F	0.8 41.7	32.9
	FRAC b	Emax FRAC D50 FRAC	0.45 18.6	-
	Tlag 1			0.401	35.1
	Tk0			1.59	31.6
	Tlag 2			22.7	-
	Ka			0.203	-
	V1			1520	-
	Q			147	89.9
	V2			2130	44.2
	CL			94.9	25.1
	σ slope (%)			20	-

a. For doses strictly inferior to 20 units 𝑭 = 𝟏, for doses superior or equal to 20 units 𝑭 = 𝟏 -𝑰𝒎𝒂𝒙 𝑭 (𝒅𝒐𝒔𝒆 -𝟐𝟎) 𝑫𝟓𝟎 𝑭 + 𝒅𝒐𝒔𝒆 -𝟐𝟎 , where dose is the administrated amount. b.

Table II .

 II Empirical estimates of Family-Wise Error Rate under H 0 for S real scenario

				FWER
	Method		C24h	C192h	AUC	CL
	Ridge regression	without correction a	13	21	15.5	9.5
	Lasso	without correction a	12	20.5	11.5	12
	HyperLasso	without correction a	10.5	17.5	11	11
	Stepwise procedure without correction a	14.5	23.5	14.5	14.5
	Ridge regression	after empirical correction b	19	21	20.5	19.5
	Lasso	after empirical correction b	18	20.5	18.5	19.5
	HyperLasso	after empirical correction b	18	17.5	18	19

Acknowledgments

Adrien Tessier received funding from Institut de Recherches Internationales Servier. The authors thank Laurent Ripoll and Bernard Walther from Institut de Recherches Internationales Servier for their advices in pharmacogenetics. The authors would also like to thank Hervé Le Nagard for the use of the computer cluster services hosted on the "Centre de Biomodélisation UMR1137".

 [31;58] Total number of true positives (TP) and false positives (FP) with their 95% confidence interval under the alternative hypothesis. On 200 simulated data sets, overall 1200 SNPs were set to impact clearance (maximum TP number).