Anne Tournadre
email: atournadre@chu-clermontferrand.fr

Bruno Pereira

Laure Gossec

Martin Soubrier

Maxime Dougados

Prof Tournadre Anne

Impact of comorbidities on fatigue in rheumatoid arthritis patients: results from a nurseled program for comorbidities management (COMEDRA)

Keywords: rheumatoid arthritis, fatigue, comorbidities

Objectives:

To analyze the factors associated with fatigue focusing on comorbidities in a large cohort of rheumatoid arthritis (RA).

Methods: Cross-sectional analyses were performed on RA patients from the French COMEDRA cohort study, a nurse-led program for comorbidities management. Fatigue was assessed using Question 3 of the Rheumatoid Arthritis Impact of Disease (RAID) score on a 0-10 numerical rating scale (NRS). Fatigue was defined as acceptable if ≤2, moderate if 3 or 4, or severe if ≥5 out of 10. Using univariate and multivariate models, the relationship between fatigue and demographics, social, disease characteristics, comorbidities (cardiovascular, infections, cancer, pulmonary, osteoporosis, and psychiatric disorders), physical activity, quality of life, and treatments was investigated.

Results: 962 patients were analyzed. The mean fatigue score was 3.8±2.7, 40% of patients reported severe fatigue. Patients had an average of 1.8 additional morbid conditions, with anxiety/depression the most common (52%). In univariate analysis, severe fatigue was more frequent in women, in patients not working, and in those with less physical activity. It was associated with disease duration and activity, mHAQ, pain, sleeping and emotional difficulties. Severe fatigue correlated with Multimorbidity index assessing the number of morbid conditions and was associated with obesity, hypertension, COPD, and anxiety/depression. In multivariate models, the risk of severe fatigue was associated with female gender, disease activity, mHAQ, current treatment with NSAIDs and biologics, multimorbidity, obesity and anxiety/depression. Conclusions: Assessment of comorbidities, psychological health and physical activity should be taken into account in order to address frequent RA-related severe fatigue.

Fatigue, along with pain, is the most frequent symptom in patients with rheumatoid arthritis (RA), affecting more than 80% of them [START_REF] Van Dartel | Association between fatigue and pain in rheumatoid arthritis: does pain precede fatigue or does fatigue precede pain?[END_REF][START_REF] Genty | Improvement of fatigue in patients with rheumatoid arthritis treated with biologics: relationship with sleep disorders, depression and clinical efficacy. A prospective, multicentre study[END_REF]. In half of patients, fatigue is severe [START_REF] Van Dartel | Association between fatigue and pain in rheumatoid arthritis: does pain precede fatigue or does fatigue precede pain?[END_REF][START_REF] Nikolaus | Fatigue and factors related to fatigue in rheumatoid arthritis: a systematic review[END_REF] and persists in almost 30-40% of cases despite conventional or biologic disease-modifying antirheumatic drugs (DMARDs) [START_REF] Genty | Improvement of fatigue in patients with rheumatoid arthritis treated with biologics: relationship with sleep disorders, depression and clinical efficacy. A prospective, multicentre study[END_REF][START_REF] Druce | Patients receiving anti-TNF therapies experience clinically important improvements in RA-related fatigue: results from the British Society for Rheumatology Biologics Register for Rheumatoid Arthritis[END_REF][START_REF] Gossec | Fatigue in rheumatoid arthritis: quantitative findings on the efficacy of tocilizumab and on factors associated with fatigue. The French multicentre prospective PEPS Study[END_REF][START_REF] Olsen | Predictors of Fatigue in Rheumatoid Arthritis Patients in Remission or in a Low Disease Activity State[END_REF]. Mechanisms are complex and multidimensional, including factors that are disease-associated (inflammation, pain, or disability), patientrelated, environment-related and medication-related [START_REF] Nikolaus | Fatigue and factors related to fatigue in rheumatoid arthritis: a systematic review[END_REF]. The majority of studies conducted to date were cross-sectional and analyzed fatigue's predictive factors. Among the main factors, associations have been demonstrated with disease activity, pain, disability, sleep, and depression, but the causal link in longitudinal studies has not been clearly established [START_REF] Nikolaus | Fatigue and factors related to fatigue in rheumatoid arthritis: a systematic review[END_REF]. Several studies have analyzed the impact of comorbidities on quality of life and disability [START_REF] Radner | Impact of comorbidity on physical function in patients with rheumatoid arthritis[END_REF][START_REF] Crilly | Relationship of EQ-5D quality of life with the presence of comorbidity and extra-articular features in patients with rheumatoid arthritis[END_REF][START_REF] Geryk | The impact of co-morbidity on health-related quality of life in rheumatoid arthritis and osteoarthritis patients[END_REF][START_REF] Mondor | The relative impact of chronic conditions and multimorbidity on health-related quality of life in Ontario long-stay home care clients[END_REF][START_REF] Van Den Hoek | Physical and Mental Functioning in Patients with Established Rheumatoid Arthritis over an 11-year Followup Period: The Role of Specific Comorbidities[END_REF], but few have focused on fatigue [START_REF] Grøn | The association of fatigue, comorbidity burden, disease activity, disability and gross domestic product in patients with rheumatoid arthritis. Results from 34 countries participating in the Quest-RA program[END_REF]. The COMorbidities, EDucation in Rheumatoid Arthritis (COMEDRA) trial is a large multicenter French cohort study providing additional data on disease-related factors, including fatigue, cognitive/emotional functioning, social aspects, comorbidities, and treatments [START_REF] Dougados | Impact of a nurseled programme on comorbidity management and impact of a patient self-assessment of disease activity on the management of rheumatoid arthritis: results of a prospective, multicentre, randomised, controlled trial (COMEDRA)[END_REF]. It offered the opportunity to explore factors associated with fatigue in a large RA cohort focusing on social aspects, comorbidities, and treatment intake.

Patients and Methods

Study population and data collection

COMEDRA is a multicenter randomized controlled trial evaluating the impact of a nurse-led program on the management of comorbidities in a French cohort of RA patients [START_REF] Dougados | Impact of a nurseled programme on comorbidity management and impact of a patient self-assessment of disease activity on the management of rheumatoid arthritis: results of a prospective, multicentre, randomised, controlled trial (COMEDRA)[END_REF]. The study comprised the recording of pre-existing comorbidities, presence of risk factors, and implementation of the recommendation for detection and management of such comorbidities. Patients with RA as defined by the 1987 American College of Rheumatology criteria, between 18 and 80 years of age, and with disease considered by the treating rheumatologist to have been stable for at least 3 months, were recruited from March 2011 to December 2011. Each patient gave his/her free, informed consent to participate, and French authorities approved the study protocol (Ethics approval: Ile-de-France III Ethics Committee, file #4-11 (B110057-30).

COMEDRA investigators collected data on socioeconomic characteristics (age, gender, employment status, and education level) and health-related risk factors (smoking status, BMI, physical activity, and alcohol consumption). Physical activity and alcohol consumption were self-reported as categorical variables 'physical activity ≥ 30 minutes by day' (yes/no) and 'alcohol consumption > 2 glasses by day' (yes/no). Furthermore, RA-related characteristics were recorded, including disease duration, tender 28-joint count, swollen 28joint count, patient and physician global assessment, ESR, CRP, presence of rheumatoid factor or anti-citrullinated protein antibodies (ACPA), erosive disease, functional impairment using the modified Health Assessment Questionnaire (mHAQ), disease impact using the Rheumatoid Arthritis Impact of Disease (RAID) score, as well as health outcome using the EQ-5D. The following data on comorbidities was recorded by the nurse: coronary heart disease, stroke, peripheral vessel disease, hypertension, diabetes, obesity, chronic kidney disease, inflammatory bowel disease, thyroid disorders, alcohol problems, depression, anxiety, fracture, bronchiectasis, chronic obstructive pulmonary disease (COPD), cancer, and history of RA-related surgeries. This was employed to calculate the Multimorbidity Index (MMI) assessing the number of morbid conditions for a patient [START_REF] Radner | Development of a multimorbidity index: Impact on quality of life using a rheumatoid arthritis cohort[END_REF]. Two types of MMI were computed, namely one enumerating comorbidities (MMI.count) and the other weighting comorbidities (MMI.weight) based on their association with health-related quality of life [START_REF] Radner | Development of a multimorbidity index: Impact on quality of life using a rheumatoid arthritis cohort[END_REF]. This index was validated on two different cohorts, and showed good performance as compared to the Charlson comorbidity index [START_REF] Radner | Development of a multimorbidity index: Impact on quality of life using a rheumatoid arthritis cohort[END_REF]. Information on treatments comprised current and past use of synthetic and biologic DMARDs, glucocorticoids, and NSAIDs. Intake of glucocorticoids was assessed using the mean prednisone dose in milligram per day during the last 3 months. Intake of NSAIDs was assessed during the last 3 months.

Assessment of fatigue

Fatigue was assessed using Question 3 of the RAID score on a 0-10 numerical rating scale (NRS) [START_REF] Gossec | Finalisation and validation of the rheumatoid arthritis impact of disease score, a patient-derived composite measure of impact of rheumatoid arthritis: a EULAR initiative[END_REF]. Respondents answered the question 'select the number that best describes how much fatigue you felt due to your rheumatoid arthritis during the last 7 days' with responses ranging from 0 (No fatigue) to 10 (Totally exhausted). Fatigue was defined as acceptable if ≤ 2, moderate if 3 or 4, or severe if ≥ 5 out of 10 (16-18).

Statistical analysis

By means of univariate analysis and multivariate models, we investigated the relationship between fatigue and demographics, social features, disease characteristics, comorbidities (cardiovascular, infections, cancer, pulmonary, osteoporosis, and psychiatric disorders), physical activity, quality of life, and treatments. All analyses were performed using Stata software (Version 13, StataCorp, College Station, TX) for a two-sided Type I error of α=5%.

Patient characteristics were expressed as mean± standard-deviation (SD) or median (interquartile range) for continuous data (assumption of normality assessed by using the Shapiro-Wilk test) and as numbers and associated percentages for categorical parameters.

The fatigue was categorized according to the predefined clinical threshold (acceptable ≤2; moderate 3-4; severe ≥ 5 out of 10) and analyzed as a categorical variable in bivariate and multivariable analyses. Because dichotomizing implies loss of information and hence loss of statistical power [START_REF] Altman | The cost of dichotomising continuous variables[END_REF], continuous variables were analyzed without categorizing them unless predefined clinical thresholds have been reported in the literature. Quantitative variables were compared between independent groups by ANOVA, or Kruskal-Wallis test if ANOVA conditions were not met (normality and homoscedasticity analyzed using the Bartlett test).

When appropriate, post-hoc tests were performed taking into account multiple comparisons (Tukey-Kramer post ANOVA and Dunn after Kruskal-Wallis). Comparisons between independent groups were carried out using the Chi-squared or Fischer's exact test for categorical variables. When appropriate, a post-hoc test was used (Marascuillo procedure).

In order to determine factors associated to fatigue considered as a 3 classes variable, a multivariable polynomial ordinal regression model was carried out using the stepwise approach (backward and forward) on covariates fixed according to univariate results (p<0.05) and clinical relevance. A particular attention has been paid to the study of multicollinearity and interactions between covariates by studying the relationships between the covariables and by evaluating the impact to add or delete variables on multivariable model. More precisely, due to this collinearity with mHAQ, pain and sleeping difficulty could not be conserved in the multivariable model including MMI count. Because mHAQ and anxiety were strongly correlated, mHAQ was not conserved in multivariable analysis including anxiety/depression but was conserved in the model with MMI count. Results were expressed as relative-risk ratios and 95% confidence intervals (95%CI) and forest plots were employed to present the results. A sensitivity analysis in patients with low disease activity (DAS28≤3.2) was performed for multivariable analyses (data not shown).

Role of the funding source

This study has been conducted thanks to a grant from the French Minister of health (PHRC 00-41).

Results

Baseline characteristics

The COMEDRA cohort included a total of 970 patients, with RAID data available for 962, which were analyzed. The mean patient age was 57.7±11.1 years, and 763 (79%) were female. The median disease duration was 11.1 years [6.2-19.1], mean DAS28 CRP 3.1±1.3 with 62 % of patients having low disease activity as defined by a DAS28 below 3.2., RAID 3±2, and mHAQ 0.40±0.45. Overall, 84% of patients were positive for RF or ACPA and 73% had erosive RA, while 673 (70%) received methotrexate, 249 (26%) NSAIDs, and 364 (38%) glucocorticoids. In total, 672 (70%) were treated with a biologic, 355 (37%) of which with TNF inhibitors, 143 (15%) with tocilizumab, 92 (9.6%) with rituximab, and 78 (8%) with abatacept.

The median MMI.count was 2 [1-3], and MMI.weight was 8 [1-10.5]. The mean MMI.count was 1.8±1.4. The most common comorbidities were anxiety/depression (52%), fracture (32%), hypertension (30%), and obesity (16%). The mean and median fatigue scores were respectively 3.8±2.7 and 4 [1; 6]. Among this stable RA population, 387 patients (40%) had severe fatigue as defined by a fatigue NRS ≥5/10 out of 10.

Relationship between fatigue, socio-demographic aspects, comorbidities, and treatment

Severe fatigue was more frequent in women, in patients not working, and in those with less physical activity (Table 1). It was associated with disease duration, disease activity, mHAQ, pain, and sleeping difficulties (Table 1). Among comorbidities, severe fatigue was associated with obesity, hypertension, COPD and anxiety/depression. In addition to data shown in table 1, no significant association was observed with cancer (n=49, p=0.15), diabetes (n=56, p=0.17) or inflammatory bowel disease (n=10, p=0.14). The Multimorbidity score (MMI count) was correlated with fatigue NRS (r=0.28, p<0.001) and patients with severe fatigue had an increased number of comorbidities (Table 1). Current treatments with NSAIDs, glucocorticoids, or biologics were associated with severe fatigue. In contrast, methotrexate use and the biologic type used did not impact the severity of fatigue.

Factors associated with fatigue level in multivariate models

The risk of being severely fatigued was associated with female gender, disease activity and mHAQ, current treatment with NSAIDs and biologics and multimorbidity (Figure 1A). Obesity and anxiety/depression were the comorbidities independently associated with severe fatigue (Figure 2A). Low physical activity, anxiety/depression and MMI count were also associated with moderate fatigue and but not female gender, obesity nor biotherapy (Figures 1,2). Results in the 2 multivariable models were similar in the subgroup of patients with low disease activity (DAS28≤3.2).

Discussion

Our results confirmed that a significant proportion of RA patients (40%) suffered from persistent severe fatigue, though the disease was considered stable and in low disease activity. Anxiety/depression, reported by more than half of patients, was strongly associated with severe fatigue. Low physical activity, others comorbidities, such as obesity, hypertension or COPD, and multimorbidity contributed to this chronic fatigue. This study confirms that RA patients are typically affected by more than one disease [START_REF] Grøn | The association of fatigue, comorbidity burden, disease activity, disability and gross domestic product in patients with rheumatoid arthritis. Results from 34 countries participating in the Quest-RA program[END_REF]. In this large RA cohort, patients had an average of 1.8 additional morbid conditions, with 50% having at least two additional conditions. MMI developed and validated in RA based on health-related quality of life [START_REF] Radner | Development of a multimorbidity index: Impact on quality of life using a rheumatoid arthritis cohort[END_REF] correlated with the intensity of fatigue. Fatigue levels were associated with female gender, disease activity, mHAQ, pain, and sleeping difficulties, in accordance with previous studies [START_REF] Van Dartel | Association between fatigue and pain in rheumatoid arthritis: does pain precede fatigue or does fatigue precede pain?[END_REF][START_REF] Nikolaus | Fatigue and factors related to fatigue in rheumatoid arthritis: a systematic review[END_REF][START_REF] Olsen | Predictors of Fatigue in Rheumatoid Arthritis Patients in Remission or in a Low Disease Activity State[END_REF][START_REF] Hifinger | In rheumatoid arthritis, country of residence has an important influence on fatigue: results from the multinational COMORA study[END_REF].

Severe fatigue was more common in patients receiving glucocorticoids, NSAIDs, and biotherapies, with similar fatigue regardless of biologic type used. Contrasting with common patient-reported side-effects [START_REF] Curtis | Patient's experience with subcutaneous and oral methotrexate for the treatment of rheumatoid arthritis[END_REF], methotrexate use was not correlated with fatigue intensity in our study. While several recent studies on biological agent use in RA reported that these drugs had a beneficial effect on fatigue [START_REF] Genty | Improvement of fatigue in patients with rheumatoid arthritis treated with biologics: relationship with sleep disorders, depression and clinical efficacy. A prospective, multicentre study[END_REF][START_REF] Druce | Patients receiving anti-TNF therapies experience clinically important improvements in RA-related fatigue: results from the British Society for Rheumatology Biologics Register for Rheumatoid Arthritis[END_REF][START_REF] Gossec | Fatigue in rheumatoid arthritis: quantitative findings on the efficacy of tocilizumab and on factors associated with fatigue. The French multicentre prospective PEPS Study[END_REF][START_REF] Chauffier | Effect of biotherapies on fatigue in rheumatoid arthritis: a systematic review of the literature and meta-analysis[END_REF], our findings indicated that current NSAIDs and biological treatments were associated with severe fatigue, which may be due to the severity of the disease even the association persisted after adjusting for disease activity. In randomized controlled trials, anti-TNF and non-anti-TNF biologic agents had similar small to moderate effects on fatigue improvement whether this was a direct or indirect effect [START_REF] Chauffier | Effect of biotherapies on fatigue in rheumatoid arthritis: a systematic review of the literature and meta-analysis[END_REF][START_REF] Almeida | Biologic interventions for fatigue in rheumatoid arthritis[END_REF]. Anti-inflammatory drugs, such as NSAIDs, glucocorticoids, and cytokine inhibitors, may directly act on pro-inflammatory cytokines involved in fatigue's biological mechanisms, in particular IL1 and IL6 [START_REF] Norheim | Biological mechanisms of chronic fatigue[END_REF]. This may impact fatigue's inflammatory mechanism, operating early in the cascade. However, other non-inflammatory pathways of fatigue may be involved at a later time point in the chronic disease process and may explain the persistent fatigue despite well-controlled disease. Disturbance of the hypothalamuspituitary-adrenal axis, cortisol production, and complex interactions with the immune system may influence fatigue [START_REF] Norheim | Biological mechanisms of chronic fatigue[END_REF]. Muscle fatigue and dysfunction during recovery induced by the production of reactive oxygen/nitrogen species due to mitochondrial dysfunction may also be involved in peripheral or physical fatigue [START_REF] Cheng | ROS/RNS and contractile function in skeletal muscle during fatigue and recovery[END_REF]. In addition, fatigue's dimensions may vary over time, and the pathogenesis of fatigue's physical and cognitive dimensions may differ, as demonstrated in cancer-related fatigue, where fatigue can persist up to 10 years after remission [START_REF] Norheim | Biological mechanisms of chronic fatigue[END_REF][START_REF] De Raaf | Inflammation and fatigue dimensions in advanced cancer patients and cancer survivors: An explorative study[END_REF].

This study displays several limitations. Because of its cross-sectional study design, the causal relationship between fatigue and associated factors could not be established. Fatigue NRS are unidimensional measures that capture only one aspect of fatigue, typically either the severity or intensity. This study can thus not analyze all dimensions of fatigue. High acceptability, completion rates, and good validity support fatigue NRS usability, although it does not allow for assessing fatigue's multidimensional aspect [START_REF] Hewlett | Measures of fatigue: Bristol Rheumatoid Arthritis Fatigue Multi-Dimensional Questionnaire (BRAF MDQ), Bristol Rheumatoid Arthritis Fatigue Numerical Rating Scales (BRAF NRS) for Severity, Effect, and Coping, Chalder Fatigue Questionnaire (CFQ), Checklist[END_REF][START_REF] Wolfe | Fatigue assessments in rheumatoid arthritis: comparative performance of visual analog scales and longer fatigue questionnaires in 7760 patients[END_REF]. There is no standardized VAS or NRS cutoff for distinguishing high and low levels of fatigue. The selected cutoff of NRS>2/10 for moderate fatigue and ≥5/10 for severe fatigue were supported by a mean VAS fatigue score in healthy controls of 20.5 ± 0.02 mm and corresponded to what is described in previous studies in RA [START_REF] Pollard | Fatigue in rheumatoid arthritis reflects pain, not disease activity[END_REF][START_REF] Hewlett | Measures of fatigue: Bristol Rheumatoid Arthritis Fatigue Multi-Dimensional Questionnaire (BRAF MDQ), Bristol Rheumatoid Arthritis Fatigue Numerical Rating Scales (BRAF NRS) for Severity, Effect, and Coping, Chalder Fatigue Questionnaire (CFQ), Checklist[END_REF][START_REF] Hifinger | In rheumatoid arthritis, country of residence has an important influence on fatigue: results from the multinational COMORA study[END_REF].

Fatigue NRS measurement was obtained using the RAID score, which was developed in an effort to combine easily, in only one questionnaire, the seven most important domains of RA impact (pain, functional capacity, fatigue, physical and emotional wellbeing, quality of sleep, and coping) without missing relevant information. The RAID questionnaire considering fatigue a core domain is very simple (seven NRS). Content and construct validity, comprehensibility, reliability, and sensitivity to change were well-documented [START_REF] Gossec | Finalisation and validation of the rheumatoid arthritis impact of disease score, a patient-derived composite measure of impact of rheumatoid arthritis: a EULAR initiative[END_REF][START_REF] Wee | Content validity of the Dutch Rheumatoid Arthritis Impact of Disease (RAID) score: results of focus group discussions in established rheumatoid arthritis patients and comparison with the International Classification of Functioning, Disability and Health core set for rheumatoid arthritis[END_REF].

This study has strengths as well involving a large cohort of patients followed up in real-life, with stable and low disease activity, and treated with various conventional and biologic DMARDs. Above all, this study allowed us to collect psychosocial and environmental variables, as well as data on comorbidities.

In conclusion, severe fatigue is a common symptom, even in patients with stable and wellcontrolled disease. In addition to anti-inflammatory medications, the possibilities of improving RA-related chronic fatigue in these patients should take into account nonpharmacological interventions, such as management of comorbidities [START_REF] Dougados | Impact of a nurseled programme on comorbidity management and impact of a patient self-assessment of disease activity on the management of rheumatoid arthritis: results of a prospective, multicentre, randomised, controlled trial (COMEDRA)[END_REF], physical activity [START_REF] Rongen-Van Dartel | Effect of Aerobic Exercise Training on Fatigue in Rheumatoid Arthritis: A Meta-Analysis[END_REF][START_REF] Durcan | The Effect of Exercise on Sleep and Fatigue in Rheumatoid Arthritis: A Randomized Controlled Study[END_REF], as well as educational programs including self-management strategies and cognitive behavioral therapy [START_REF] Hewlett | Self-management of fatigue in rheumatoid arthritis: a randomised controlled trial of group cognitive-behavioural therapy[END_REF].

Acknowledgment and funding.

Table1. Relationships between fatigue and socio-demographic aspects, comorbidities, and treatments. Table1. Relationships between fatigue and socio-demographic aspects, comorbidities, and treatments.

Figure 1 .

 1 Figure 1. Forest plot showing the overall Risk ratio in multiple regression model including multimorbidity score (MMI.count) for severe fatigue (A) and moderate fatigue (B).

Figure 2 .

 2 Figure 2. Forest plot showing the Risk ratio in multiple regression model including each comorbidity for severe (A) or moderate fatigue (B).

We thank the patients and investigators who participated in the study : AP-HP (Cochin, Lariboisiére, Pitié-Salpêtrière), CHU Bordeaux, CHU Brest, CHU Clermont-Ferrand, CHU Grenoble, CH Le Mans, CHU Lille, CHU Marseille, CHU Montpellier, CHU Nancy, CHU Rennes, CHU Rouen, CHU Strasbourg, CHU Toulouse.

Competing interests: None.

Patient consent Obtained.

Ethics approval: Ile-de-France III Ethics Committee, file #4-11 (B110057-30).

Trial registration number Clinical trials NCT #01315652


Bertrand Julie
email: <julie.bertrand@inserm.fr>

Comets Emmanuelle

Julie Bertrand

Comparison of model-based tests and selection strategies to detect genetic polymorphisms influencing pharmacokinetic parameters

Keywords: Adult, Algorithms, Analysis of Variance, Bayes Theorem, Computer Simulation, Drug Monitoring, Female, HIV Infections, drug therapy, genetics, metabolism, HIV Protease Inhibitors, pharmacokinetics, therapeutic use, Humans, Indinavir, pharmacokinetics, therapeutic use, Likelihood Functions, Male, Models, Biological, Models, Statistical, Nonlinear Dynamics, P-Glycoprotein, genetics, metabolism, Pilot Projects, Polymorphism, Single Nucleotide, Software Pharmacokinetics, Genetics, nonlinear mixed effects model, test, model selection, NONMEM

We evaluate by simulation three model-based methods to test the influence of a single nucleotide polymorphism on a pharmacokinetic parameter of a drug: ANOVA on the empirical Bayes estimates of the individual parameters, likelihood ratio test between models with and without genetic covariate, and Wald tests on the parameters of the model with covariate. Analyses are performed using the FO and FOCE method implemented in the NONMEM software. We compare several approaches for model selection based on tests and global criteria. We illustrate the results with pharmacokinetic data on indinavir from HIV-positive patients included in COPHAR 2-ANRS 111 to study the gene effect prospectively. Only the tests based on the EBE obtain an empirical type I error close to the expected 5 . The approximation made with the FO algorithm results in a significant inflation of the type I error of the LRT % and Wald tests.

INTRODUCTION

Pharmacokinetics studies the time course of a drug in the body (

). The variability in the [START_REF] Gabrielsson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] pharmacokinetics of drugs when administered to different subjects is often important and should be studied to improve the use of drugs, avoid toxic events and allow individualization of therapy. The contribution of genetic factors to this variability is potentially important ( ). Some genes have already been the subject of much attention, for example the ABCB1 gene coding for the P [START_REF] Licinio | Pharmacogenomics The search for individualized therapy[END_REF] glycoprotein (P-gP) found on the main exchange barriers (

). The involvement of this protein in drug absorption [START_REF] Marzolini | Polymorphisms in Human MDR1 (P-glycoprotein): Recent advances and clinical relevance[END_REF] processes has been demonstrated directly in animals and indirectly in humans. For example, co-administration of paclitaxel and cyclosporine, respectively a substrate and an inhibitor of P-gP, has been shown to increase the bioavailability of paclitaxel (Meerum et al., ). More recently, Yamaguchi and colleagues demonstrated the impact of ABCB1 polymorphism on paclitaxel pharmacokinetics in 1999 patients with ovarian cancer (

). The ABCB1 gene is composed of 209 kb ( ), and to date 28 [START_REF] Yamaguchi | Genetic variation in ABCB1 influences paclitaxel pharmacokinetics in Japanese patients with ovarian cancer[END_REF][START_REF] Bodor | Characterization of the human MDR1 gene[END_REF] polymorphisms concerning a nucleotide modification have been described, including several with potential clinical impact (Sakaeda et al., ). In the present study, we consider data from patients treated with indinavir in the COPHAR2-ANRS11 study investigating the 2002 benefit of early therapeutic drug monitoring in anti-retroviral therapy. Since the ABCB1 genetic polymorphism was found to have an influence on the pharmacokinetics of protease inhibitors ( ), patients were genotyped for two exons (exon 21 and 26) of [START_REF] Fellay | Response to antiretroviral treatment in HIV-1-infected individuals with allelic variants of the multidrug resistance transporter gene MDR1: a pharmacogenetic study[END_REF] the ABCB1 gene.

The influence of genetic polymorphism on concentration data is usually analysed using non-compartmental analysis (NCA). The area under the concentration time curve as well as model-independent parameters are calculated for each individual concentration profile versus using, for example, the log-trapezoidal method (

) and individual pharmacokinetic (or PK) parameters are [START_REF] Gabrielsson | Pharmacokinetic and Pharmacodynamic Data Analysis: Concepts and Applications[END_REF] then compared between the different genotype groups using ANOVA (

). More sophisticated approaches using [START_REF] Inomata | CYP2C19 genotype affects diazepam pharmacokinetics and emergence from general anesthesia[END_REF] nonlinear PK models in individual regression ( ) or mixed effects models ( ) have also been applied to [START_REF] Min | CYP3A5 polymorphism and the ethnic differences in cyclosporine pharmacokinetics in healthy subjects[END_REF][START_REF] Taguchi | Nonlinear mixed effects model analysis of the pharmacokinetics of routinely administered bepridil in japanese patients with arrhythmias[END_REF] genetic data in PK studies. Nonlinear mixed effects models (NLMEM) require fewer blood samples in each patient and can be important in special populations such as patients with acute disease or neonates, for whom extensive sampling is obviously impractical.

Various methods can be used to include pharmacogenetic information in NLMEM ( ). Empirical Bayes estimates Comets et al., 2007 (EBE) of the individual parameters can be computed from the fit of a model with no covariate and compared between the different genotype groups with an ANOVA to test the genetic polymorphism effect (

). Another approach is to perform [START_REF] Henningsson | Association of CYP2C8, CYP3A4, CYP3A5, and ABCB1 polymorphisms with the pharmacokinetics of paclitaxel[END_REF] Wald tests on the estimates of the gene effect coefficients in the covariate model ( ), while stepwise model building is [START_REF] Kerbusch | Population pharmacokinetic modelling of darifenacin and its hydroxylated metabolite using pooled data, incorporating saturable first pass metabolism, CYP2D6 genotype and formulation-dependant bio availability[END_REF] frequently based on the comparison of models with and without gene effect using the likelihood ratio test (LRT) (

). [START_REF] Mamiya | CYP2C9 polymorphism effect on phenobarbiton kinetics in Japanese patients with epilepsy: analysis by population phamacokinetics[END_REF] The first purpose of this work is to assess the statistical properties of these three different strategies to test for a gene effect through a simulation study. The setting for the simulation study is based on the COPHAR2-ANRS11 study. We have also taken the impact of estimation methods into account, comparing results from the FO and FOCE methods implemented in the NONMEM (nonlinear mixed effect model) software (

). The type I error of the tests is evaluated using simulations under the null hypothesis H of [START_REF] Sheiner | NONMEM Version 5.1. San Francisco University of California[END_REF] no gene effect for two designs (40 patients as in the study, 200 patients to examine the influence of sample size), while the power is compared using simulations under an alternative hypothesis H .

1

In the same simulation study, we also examine model selection strategies, in which the aim is to choose the best covariate model for the gene effect. We compare test-based strategies with general selection criteria that have been proposed, like the Bayes information criterion (BIC) (

) and the Akaike information criterion (AIC) ( ). According to the parsimony principle, these [START_REF] Schwartz | Estimating the dimension of a model[END_REF][START_REF] Akaike | A new look at the statistical model identification[END_REF] criteria penalize the log-likelihood by the number of model parameters which limits overfitting. Alternative criteria derived from the AIC have been developed in order to deal with small samples or to improve AIC asymptotic consistency. The properties of these criteria have been studied in generalized linear models ( ), like the trend to select an over-parameterized model with AIC [START_REF] Burnham | Model Selection and Multimodel. A practical Information-Theoric approach[END_REF] and under-parameterized model with BIC. AIC and BIC have been widely used in population PK and PK/PD studies during the past decade. In the present paper, we compare through a simulation study their properties for nonlinear mixed effects model selection.

In this paper, we first introduce the model, the notations, the different tests and model building strategies under study and the estimation methods. Then we describe the case study and the simulation study based on its design. Next we present the results of the simulation study and the application to the real data. Finally we discuss our findings.

METHODS

Model and notations

Let the function denote the PK model, which depends nonlinearly on its parameters . The concentration at time for subject f θ y i , j t i , j i 1, , and measurement 1, , is given by:

= … N j = … n i
where is the vector of PK parameters for the i individual and the residual error, assumed to be normally distributed with zero

θ i P th ε i , j
mean and variance . Here we assume a proportional error model:

σ i , j 2
We assume a multivariate log-normal distribution for the vector of individual parameters, are then expressed as:

θ i
with the vector of fixed effects and the vector of random effects. The random effects are assumed to be independent of and

μ b i b i ε i , j
normally distributed with zero mean and variance matrix . In this work, we used a diagonal with interindividual variance of the p Ω Ω th parameter , however an unspecified positive definite matrix can be assumed.

G i β G i = β 0 β 1 β 2 G i =
TT, respectively. We write the model for the genetic polymorphism effect in subject as: i

We assume that CC is the reference class so that 1.

β 0 =
Let model M be the model in the absence of genetic polymorphism effect: 1 , (CC CT TT). There are three base

{β 0 = β 1 = β 2 = } = =
models including the gene effect as a covariate: a complete model where the mean of the parameter is different in the three groups, M : mult {

1, 1 , (CC CT TT) and two reduced models, M : 1, 1 , (CC CT TT) and M : 1,

β 0 = β 1 ≠ β 2 ≠ } ≠ ≠ recessive {β 0 = β 1 = β 2 ≠ } = ≠ dominant {β 0 = β 1 = β 2 1 , (CC CT TT). ≠ } ≠ =
In the following, corresponds to the model log-likelihood, represents the number of population model parameters (mean PK L P pop parameters, covariate coefficients, variances and error model parameters), the sample size and the total number of observations. The N n tot standard error and covariance of the parameter estimates are abbreviated by SE and cov, respectively.

Tests for genetic polymorphism effect

In this section, we examine different tests based on NLMEM that can be used to test the existence of a genetic polymorphism effect on one parameter of a PK population model.

Analysis of variance (ANOVA)

Data are analysed with a model with no covariates (M ) and Empirical Bayes estimates (EBE) of the individual PK parameters are base computed. A one-way analysis of variance is used to detect differences in these EBE between the different genetic groups.

Wald test

The data are analysed using the complete model M , and the significance of the parameters is assessed using a global Wald test: mult where represents the variance covariance matrix of parameters and . The statistic is compared to the critical value of a

Σ β 1 β 2 W χ 2
with two degrees of freedom.

Likelihood Ratio Test (LRT)

The third test relies on the comparison between the model with no covariate effect M and the complete model M . These two number of population parameters between the two models.

Strategies for model building

In this section, we examine a second aspect, which is model selection. The different strategies under study provide decision rules for covariate inclusion in order to get the best model. represents in a diagram the decision path of the three strategies using tests. Figure 1 When performing multiple tests, a correction which preserves the false discovery rate is used ( ). [START_REF] Benjamini | Controlling the false discovery rate: a practical and powerful approach to multiple testing[END_REF]

Selection based on EBE

If the test that compares means of EBE obtained with M between CC, CT and TT using an ANOVA is significant, the means of the base EBE obtained with M are compared between the CC and the CT, on one hand, and between the CT and the TT, on the other hand, using dominant

Selection based on information criteria

As an alternative, criteria such as AIC or BIC can be used for model selection. Here, the model with the lowest criterion is chosen.

Both of these criteria balance the log-likelihood with the number of parameters in the population model. According the parsimony principle, a simpler model is preferred for equivalent information gain. The AIC is written as:

The expression of the BIC involves the total number of observations:

The BIC is related to the Bayes factor which is a measure of the strength of evidence in favour of a given model and can thus be used to quantify model uncertainty.

Other criteria have been defined, derived from those above: the corrected AIC (AICc) defined by Akaike for small samples i.e. when n / < 40 ( ) and the consistent AIC (CAIC) ( ). They are derived from AIC and involve the total number tot P pop [START_REF] Sugiura | Further analysis of the data by Akaike s information criterion and the finite corrections[END_REF][START_REF] Bozdogan | Model selection and Akaike s information criteria (AIC): the general theory and its analytical extensions[END_REF] of observations:

Another formulation for BIC has also been proposed, BICc, where is replaced by the number of subjects ( ): n tot Raftery, 1995

Estimation methods

In NLMEM, the parameters and their standard error are mainly estimated using maximum likelihood. However the likelihood function for these models is expressed as an integral and has no analytical solution. Specific algorithms have therefore been proposed to perform the maximization. The approach most frequently used relies on first-order approximations of the likelihood function. This approach has been implemented in the NONMEM software version V which is the most frequently used software in PK/PD analyses. In NONMEM two approximations can be used: the first-order method (FO), relying on a first-order linearization of the likelihood function around 0, and b i = the first-order conditional estimation method (FOCE), relying on a first-order linearization of the likelihood function around the estimates of the individual random effects. Although the FO method suffers from bias and lack of precision it is still frequently used as it is faster and has less convergence problems than FOCE. In this paper, we compare both estimation methods. More precisely, we use the FOCE with interaction method, allowing possible interaction between inter-and intra-variability.

Real data and simulation study

Real data

We illustrate the different approaches on data from a PK sub-study of the COPHAR 2-ANRS 111 study, a multicentre non-comparative pilot trial of early therapeutic drug monitoring in HIV positive patients na ve of treatment. The objective of the trial was ï to assess the benefit of a pharmacological intervention after measurement of trough plasma concentrations of protease inhibitors (Mentr et é ). We focus on the PK sub-study from the group of patients receiving indinavir boosted with ritonavir. Patients were genotyped al., 2005 for the ABCB1 exons 21 and 26 to investigate genetic polymorphism impact on pharmacokinetics of the protease inhibitors, which are well-known substrates of the P-gp (

). [START_REF] Fellay | Response to antiretroviral treatment in HIV-1-infected individuals with allelic variants of the multidrug resistance transporter gene MDR1: a pharmacogenetic study[END_REF] Forty-two patients were included, one patient withdrew from the study and one switched to another protease inhibitor during the first week of treatment. We therefore obtained PK data for 40 patients (27 men, 13 women) with an average age of 36.5 years. PK profiles were determined at 1, 3, 6 and 12 h after administration of the drug at a date two weeks after the treatment onset. One patient had missing information for the two genotypes, and ABCB1 exon 26 genotype was missing in two other patients.

Simulated data

The design used in the simulation mimics that of the application data set. We simulate PK studies of N 40 (equivalent to the where represents the bioavailability, the absorption rate, the elimination rate, the volume of distribution and 12 h the time

F k a k V τ =
between two doses. We use the following population parameters: / 102 L with interindividual variability 41.3 ,

1.4 h V F = ω V = % k a = 1 - with
113 , 0.2 h with 26.4 obtained with a preliminary analysis of the real data using the FO method in NONMEM ω k a

= % k = 1 - ω k = %
and a residual error of 20 . The measurements below the quantification limit (BQL) are treated in the analysis using a standard approach:

% the first value in a series of BQL was set to LOQ/2 and the remaining values were discarded ( ). [START_REF] Beal | Conditioning on certain random events associated with statistical variability in PK/PD[END_REF] Under H , 1000 data sets are simulated with a design of N 40 patients and 1000 with a design of N 200 patients. We simulate a ) is , or if GG, GT or TT. Under H we set (G residual error, generated from a normal distribution (0, ( , )), to each predicted concentration to obtain the simulated N σ 2 f θ t i , j concentration.

β 26i β 0 β 1 β 2 = δ G 21i δ 0 δ 1 δ 2 G 21i = 1 β ) 1, 1.2, 1.
Simulations were performed using the statistical software R (R ) running under Linux (Red Hat 9.0).

Development Core Team 2006

Evaluation of the tests for genetic polymorphism effect

In the first step, each of the three tests presented above is applied to detect the effect of the exon 26 polymorphism on the 1000 data sets simulated under H for each design (40 and 200 patients). Tests are performed on estimations obtained for each data set using FO and 0 FOCE in NONMEM. The type I error for each analysis is defined as the percentage of data sets where the corresponding test was significant. The expected prediction interval with 1000 simulations and a value of 5 is 3.7; 6.3 . To ensure a type I error of 5 , we

% [ ] %
define a correction threshold as the 5 percentile of the distribution of the p-values of the test under H .

th 0

In the second step, for the design with 40 patients, the same tests are performed using the 1000 data sets simulated under H . The 1 power is defined as the percentage of data sets where the corresponding test was significant. We use the corrected threshold to compute the corrected powers. This allows comparison of the different tests even if the type I error is different from 5 .

% Data sets with a group defined by the genotype for exon 26 with less than two patients are discarded from the analysis. It should be emphasized that the number of simulated data sets for which the tests can be applied are often less than 1000. Indeed, the algorithms FO and FOCE used in the NONMEM software are sometimes unable to converge on a data set. In that case, estimates are not available and no tests can be applied. Similarly, even when the linearization algorithms achieve convergence, the variance covariance matrix is not always available. More precisely, to perform the ANOVA on the EBE only the convergence of model M is required with or without the SE base estimates. The Wald test requires the parameter estimation error to compute the test statistics, thus not only the convergence of model M is required but the variance covariance matrix must also be available. The convergence of both M and M is necessary to apply mult base mult the LRT but the covariance step is not needed to succeed. The three tests were also evaluated on a subset of data sets fulfilling all the conditions listed in this paragraph using FO or FOCE.

Evaluation of the strategies for model building

The As for the evaluation of tests, the analyses are not always performed on the 1000 data sets. The conditions to perform model building following an ANOVA are identical to those described for the test. To choose the best model using Wald tests, the convergence and the covariance step are required for the three models including the gene effect. For the selection with the LRT, convergence of the four models is required and for the criteria based selection, at least one model must have achieved convergence. The different strategies were also evaluated on samples satisfying all conditions with FO and FOCE.

Application to real data

The indinavir concentrations are analysed with the same PK model as in the simulation study. The estimation method is chosen based on performances in the simulation study. Both the exons 26 and 21 of the ABCB1 gene are investigated using the three tests and the eight model building strategies.

RESULTS

Type I error and power of the detection tests

For the design with N 40, three simulated data sets (one under H and 2 under H ) are discarded from the analysis due to a group of = 0 1 less than two patients with the TT genotype. Also, FOCE encounters many more convergence problems than FO; , the number of data K sets among the 1000 on which the test is performed, is always lower with FOCE.

Type I error estimates of the three tests performed are shown in . For the design with N 40, the ANOVA type I errors do not Table 1 = significantly differ from 5 with both estimation algorithms. The LRT type I error estimates show a slight significant increase with FOCE % and increase by ten times using FO. For the Wald test, there is a rather important significant increase for FOCE and again inflation by four times with FO. For the design with 200 patients, the LRT attain a type I error non-significantly different from 5 and the Wald test still % has an estimate slightly superior to the nominal level with FOCE. The large increase of the type I error remains the same with FO for the LRT and the Wald test.

Estimates are given of the power of the tests using both estimation algorithms for the ANOVA, the Wald test and the LRT, for the design with 40 patients, in . For each test and each estimation method, the corrected power is computed using the corresponding Table 2 empirical threshold in order to maintain a type I error of 5 . The corrected power for LRT and Wald tests based on the FO method were % low, but the inflated type I error in these situations already show that this estimation method is poor. The powers for the ANOVA (using FO or FOCE) and the LRT for FOCE are around 70 , but the power of the Wald approach for FOCE is much lower (25 ). We explored % %

the FOCE outputs to understand this last result. In fact, with the FOCE algorithm in NONMEM, we observe correlations between the estimates of the gene effect coefficients and their estimation errors. This relationship leads to decreased values of the Wald statistic and therefore reduces the power to detect a genetic polymorphism effect. The same results were obtained considering only the subset of data sets fulfilling all convergence conditions for both FO and FOCE.

Model selection strategies

The results of the model selection strategies applied to the data sets simulated under H with the design of 40 patients are presented in 0 . Again, the performances of the FO algorithm are unsatisfactory except for the strategies based on EBE. With FOCE, both model Table 3 selection strategies based on EBE and LRT select the correct model M in about 90 of instances whereas the Wald approach selects a base % model with a gene effect in around 15 of the data sets. Using selection criteria, the AIC and AICc obtain the worst performances, % selecting a model with a gene effect in 57.8 and 51.7 of the data sets respectively. BICc shows performance close to that of the % % approach based on Wald tests, while CAIC (4.7 ) and BIC (6.7 ) select a model with a gene effect in about 5 of the data sets.

% % %

The results on data sets with N 40 simulated under H are given in . The results of the strategies using the LRT, Wald tests or

= 1
Table 4 a criterion with the FO estimation method are not presented because of the poor properties under H . A model with a gene effect is 0 selected for about 70 of the data sets using CAIC, for about 80 of the data sets using the different model-building strategies based on % % tests or BIC, for about 90 of the data sets using BICc, and for about 99 of the data sets using AIC and AICc. However, for the latter, selecting a model with a gene effect, it is not always the correct one. The performances are similar using the sample of data sets for which all conditions were put together with FO and FOCE.

represents a summary of these results, placing each strategy on a bi-dimensional plan; the ability to select a model with a Figure 2 gene effect under H versus the trend not to select M under H . The CAIC is the test which has the best properties under H , but the 1 base 0 0

BIC more often selects a model with covariate under H while it is only slightly less conservative under H . We also note that the LRT, the 1 0 BIC and the ANOVA are clustered together, and thus offer similar compromises.

Application

represents a spaghetti plot of concentrations versus time for the 37 patients from the indinavir arm of the COPHAR2 study, Figure 3 sorted by genotype classes for exon 26. Concentrations show an important interindividual variability and only three patients were mutant homozygotes for this polymorphism.

As the simulation study has shown poor performances with the FO algorithm, we estimated the model parameters using the FOCE with interaction algorithm. Therefore, the estimates are different from the simulated values obtained in the preliminary analysis. In addition, the estimate of the interindividual variability for k was very small and we fixed it to 0 (no variability). The model with no covariate had an absorption constant of 0.8 h with an important interindividual variability of 70.3 , an elimination constant of 0.

DISCUSSION

In this work, we evaluate several statistical tests and model selection strategies using nonlinear mixed effects models to analyse the impact of a genetic polymorphism on one PK parameter through simulation. We also study the impact of the estimation algorithms in settings. Of course if the study was designed specifically to detect a gene effect the sample size could be increased to ensure a higher (more ethical) power. Our objective here was to compare different methods and we used the data from the COPHAR2-ANRS 111 trial to provide the settings for the PK simulations, as we have analysed these data thereafter. The simulated effect of gene was chosen to be consistent with the literature, though in the COPHAR2 trial the genotype distribution is slightly unbalanced. The reduced power of the either for numerical reasons or because they had to be terminated and we could not obtain an estimate of the variance-covariance matrix in 32 (3.2 ) and 38 (3.8 ) runs, respectively. This could have been improved partly by the use of different initial conditions, nonetheless the % % results were identical using the sample of data sets fulfilling convergence conditions for the three tests using both estimation methods.

Among the different implementations of both algorithms available, we chose to use NONMEM as it is the most popular tool in the pharmaceutical industry. Our results can be extended to the FO method implemented in SAS as it computes the same likelihood function up to a constant ( ), although it is not the same algorithm ( ). Similarly, FOCE with interaction and the nlme method [START_REF] Wang | Derivation of various NONMEM estimation methods[END_REF][START_REF] Roe | Comparison of population pharmacokinetic modeling methods using simulated data: results from the Population Modeling Workgroup[END_REF] with the varConstPower option implemented in Splus should give the same results, save for the lower accuracy of the Splus approximation ( , ). The simulations and analyses in the present paper have been performed using version 5.1 of Girard and Mentr , 2005 é Wang, 2007 NONMEM. A new version of NONMEM, version 6, has been released in December 2006, and was implemented in our department after the major part of the analyses had already been run, so that we kept NONMEM 5.1 for this study. The FOCE routine in NONMEM 6 has been rewritten and should provide more stable runs, which may reduce the convergence issues we have found in the present study. However, in our preliminary results with using NONMEM 6 on a subset of the simulated datasets, an increased number of runs failed the covariance-step. Rather than reporting partial results in the present study, we will investigate this matter in a subsequent work. . Indeed, in our simulation conditions, M and M close and the power to detect a difference between these two models is recessive recessive mult much lower than the power to detect a difference between M and M . As a side-note, in the Bayesian literature, it is usual also to base mult consider not only the best model (i.e. the model with the lowest criterion), but also models close to the best model ( ). [START_REF] Raftery | Bayesian Model Selection in Social Research (with discussion)[END_REF] However in our study, the simulation model is rarely close in this sense to the best model.

Finally, in our study, for the design with 40 patients, the ANOVA on the EBE is the only test that maintained a 5 type I error, as well % as a good power. It should be noted that we simulated a sufficient number of samples per patient with respect to the number of model parameters. In sparse sample situations, regression to the mean is known to occur with EBE ( ), which could Panhard and Mentr , 2005 é result in lower power to detect differences between genotype groups. We plan to test this hypothesis in a subsequent study with a sparse design including only two samples per patient.

We also confirm that both AIC and AICc should not be used for model building but further studies are required to provide recommendations on the other selection strategies. Regarding the estimation methods, if FO can still be used in covariate screening on the EBE, one should avoid performing model building with the LRT or Wald tests based on results from this algorithm.

We illustrate the different approaches using data from the indinavir PK sub-study of the COPHAR2 ANRS-111 trial. The PK model has already been described but no gene effect has been investigated using a population approach to date. The estimated parameters are in accordance with estimations obtained in other studies ( , , ). The Wald test is [START_REF] Csajka | Population pharmacokinetics of indinavir in patients infected with human immunodeficiency virus[END_REF]Goujard et al., 2005[START_REF] Kappelhoff | Population pharmacokinetics of indinavir alone and in combination with ritonavir in HIV-1-infected patients[END_REF] the only test to detect an influence of ABCB1 exon 26 on the volume of distribution. Considering the corrected threshold provided by the simulation study, we can probably ascribe this discrepant result to the inflated type I error. Similarly, it is the three strategies with the highest percentage of data sets where M is not selected under H that select a model with the ABCB1 exon 26 polymorphism as base 0 covariate. The polymorphism on the exon 26 of the ABCB1 gene has been shown to impact on plasma concentrations of nelfinavir, another protease inhibitor ( ). However, this work agrees with another study ( ), where no effect of [START_REF] Fellay | Response to antiretroviral treatment in HIV-1-infected individuals with allelic variants of the multidrug resistance transporter gene MDR1: a pharmacogenetic study[END_REF][START_REF] Verstuyft | Absence of association between MDR1 genetic polymorphisms, indinavir pharmacokinetics and response to HAART[END_REF] the polymorphisms from ABCB1 exon 26 and 21 on indinavir bioavailability was found.

Another extension of this work would be to simulate under H various levels of the gene effects, which would provide more 1 information on the relationship between the strength of the genetic polymorphism effect and the power. Designing an optimal sampling schedule for testing a gene effect with a given power is also an interesting challenge and the extension of the PFIM software for design optimization in the case of models with covariates could be used ( , ). Moreover in the genetic Retout andMentre, 2003 Retout et al., 2007 framework one should keep in mind the complex pathway leading from DNA to metabolic activity, which is usually controlled by more than one exon. Another perspective would therefore be to analyse the influence of the haplotypes (

), since such a [START_REF] Innocentia | Haplotypes of variants in the UDP-glucuronosyltransferase 1A9 and 1A1 genes[END_REF] classification seems to be more relevant at the DNA level.

Pharmacogenetic studies using NLMEM have many advantages as fewer samples are required to estimate parameters with a biological meaning. The current literature presents a wide array of methods for covariate selection using NLMEM. We show in this study that methods using EBE are not only efficient in data exploration but also in model selection on data sets with enough samples per patient. We also emphasize that using estimation algorithms based on likelihood linearization, LRT type I error is inflated, thus one has to perform simulations or work with large data sets. Finally, the problems of the FO and FOCE algorithm in terms of convergence and bias are an incentive to use more recent estimation methods ( ). [START_REF] Samson | The SAEM algorithm for group comparison test in longitudinal data analysis based in nonlinear mixed-effects model[END_REF] Ackowledgements:

We would like to thank the COPHAR 2-ANRS 111 scientific committee (investigators: Pr. D. Salmon and Dr X. Duval, pharmacology: Pr JM. Tr luyer, methodology: Pr F. Mentr ) for giving us access to the PK data of the indinavir arm in order to build our simulations and to é é illustrate our topic. We would also like to thank the IFR02 of INSERM and Herv Le Nagard for the use of the centre de biomod lisation .

é " é " Results obtained with FO are not presented for these strategies because of their poor performance under H0 ( ). table 3

  let us assume that our aim is to detect the effect of a single nucleotide polymorphism (SNP) on one PK parameter, for instance the p component of . Let C denote the wild-type allele and T the mutant. An individual can be one of three genotypes (eg. th ; θ θ p CC for the wild-type homozygotes, CT or TT). Let denote the genotype for subject i, and let ( ) , or for CC, CT or

  compared with the LRT. The test statistic 2() is compared with a with two degrees of freedom S LR = -L base -

  base t-tests. If none or one test is significant, the one with the lower p-value leads to model selection, or else M is selected. For example, if mult the test comparing the CC to the CT has the lower p-value, M is chosen. dominant If the global test is non-significant, the means of the EBE obtained with M are compared using t-tests opposing: i) the CC to the CT base and TT put together, and ii) the TT to the CT and CC put together. If one or both tests are significant, the model corresponding to the test with the lower p-value is selected, else M is selected. For example if the test comparing the CC to the group formed by the CT and the base TT has the lower p-value, M is chosen. dominant Selection based on Wald test If the global Wald test comparing to a with two degrees of freedom is significant, two Wald tests on coefficients estimated in M W χ 2 are then realized; comparing both and with a with one degree of freedom. If none or one test is significant, the one with mult χ 2 the lower p-value leads to model selection, or else M is kept. As an example, if has the lower p-value, M is chosen. mult dominantIf the global test is non-significant, data are analysed using the two models with the gene effect in two classes and Wald tests are performed on the coefficients estimated in If one or both tests are significant, the model is chosen based on the test with the lower p-value, or else M is conservednone or one test is significant, the test with the higher p-value leads the model choice, or else M is selected. As an recessive mult example, if the test comparing M with M obtains the higher p-value, M is selected. mult dominant dominant If the global test is non-significant, M is compared with M and M using LRT. Then, if one or both tests are significant, base dominant recessive the model is chosen based on the test with the lower p-value, or else M is selected. For example, if the test comparing M to M base base dominant obtains the lower p-value, M is selected.

=

  application sample size) and N 200 patients with four samples (1, 3, 6, 12 h after dose) at steady state. The same bid doses of 400 mg for = indinavir and 100 mg for ritonavir are assumed for all patients. The concentrations are simulated using the steady-state one-compartment model with first-order absorption and elimination that was used to model the indinavir concentrations in the COPHAR 1-ANRS 102 study():[START_REF] Brendel | Population pharmacokinetic analysis of indinavir in HIV-infected patient treated with a stable antiretroviral therapy[END_REF]

  on two exons located on the same chromosome. The bioavailability is assumed to depend on the diplotype. Because F the model is parameterized as / , and , this is equivalent to assuming that / depends on the diplotype. The distribution of the exons V F k a k V F mimics that of exon 26 and exon 21 of the ABCB1 gene as reported by Sakaeda and colleagues ( ) we note C and G Sakaeda et al., 2002respectively the wild-type allele for the 2 exons and T the mutant allele. With those properties, for data sets with 40 patients we expect on average 9 individuals with a CC genotype, 18 with CT genotype and 13 with TT genotype for exon 26. Under the alternative hypothesis, we assume the following effect of the two polymorphisms (modified from ), where equation (4) G 26i denotes the genotype for the exon 26 and the genotype for the exon 21: G 21i where (G ) is , or if CC, CT or TT as previously and (

  6 and ( ) 1, 1.1, 1.3. These values were chosen to provide a good power for the detection of one SNP effect in the 26 = δ G 21 = context of the simulations while remaining consistent with results found in the literature concerning the effect of ABCB1 polymorphism. With those genetic coefficient effects and the distribution from the literature, we simulated 100 data sets with the N 40 design. When we = computed the ( / ) EBE from M and performed bilateral t-tests to compare the wild homozygotes and the mutant homozygotes for V F i base the exon 26, we obtained a power of 80 (). %[START_REF] Machin | Sample size tables for clinical studies[END_REF] The simulation of each data set is performed as follows. The set of possible genotypes for exon 26 and exon 21 is S CC-GG, CC-GT, ={ CC-TT, CT-GG, CT-GT, CT-TT, TT-GG, TT-GT, TT-TT with corresponding simulated frequencies f 0.2, 0.02, 0.02, 0.05, 0.38, 0.05, } ={ 0.04, 0.04, 0.2 . For each individual the genotypes are drawn from this distribution. Under H , both genotypes condition the value of the } 1 fixed effects for V/F according to . Then we simulate a random effect vector from a normal distribution (0

-

  of distribution of 99.3 L with an interindividual variability of 47 . All the estimation errors were below 20 for the fixed % % effects and below 40 for the variances. % For the influence of the ABCB1 exon 26 on the indinavir volume of distribution, the ANOVA and the LRT were non-significant (P = 0.7 and 0.2 respectively). The global Wald test obtained a p-value of 0.02, however the corrected threshold defined in the simulation study for this test is 7.5.10 . The model with no covariate was chosen using the selection strategies based on the EBE, or on the LRT, using 3 CAIC, BIC and BICc, whereas the strategy based on Wald tests, AIC and AICc selected M . No influence of the ABCB1 exon 21 on recessive indinavir volume of distribution was detected using any of the tests or selection strategies under study.

é

  NONMEM, comparing the two first-order approximations most widely used FO and FOCE, which linearize the model function respectively around the random effects equal to 0 and around the individual estimates of the random effects. Although the FO method has been shown to suffer from various problems, it is still used because the FOCE algorithm is known for numerical difficulties and for its slowness. Using the FO algorithm we observe unsatisfactory performances for all the tests and model selection strategies with the exception of methods based on the EBE. The linearization of the likelihood function around the fixed effects leads to type I error inflation With FOCE, there is a significant increase in the type I error of the LRT and the Wald approach with a design including 40 patients. This increase has already been described for the LRT ( , ) and has been also shown for theComets and Mentr , 2001 é W hlby et al., 2001 ä Wald test ( ). The design with 200 patients is closer to asymptotic conditions and shows as expected a correct Panhard and Mentr , 2005 é type I error. Performing simulations under H can be used to correct the threshold for the test under H as we did in the present study. We 0 1 observe a power around 70 for the tests using ANOVA and LRT which is close to the power of 80 expected from the simulation % %

  Wald approach and the unsatisfactory efficacy of the model selection strategy based on the Wald tests could result partly from a wrong estimation of the standard errors due to the log-likelihood function linearization. Indeed, we observe with FOCE that in the simulations the estimation errors are highly correlated to their estimates. Because the Wald statistic is based on the ratio of the estimates to the estimation error, this could explain the poor performance of the Wald tests under H . Finally, FOCE met with convergence problems: on M and M 1 0 which involved no mathematical complexities for N 40 under H , 35 (3.5 ) and 48 (4.8 ) runs did not achieve convergence,

  Model selection strategies based on tests have, to a certain extent, a high rate of false inclusion under H (over 10 ) which could 0 % result from the unconnected multiple model comparisons. Further there is no simple way to correct the model selection under H by taking 1 into account the simulations under H . Consequently, we have to be cautious about assessing the performance under H selection based on criteria, except for the performances of the BIC and the BICc not showing any trend to conservatism, under H our results agree with the literature. The acceptable performance of the CAIC has to be noted, as well as the very 0 poor performance of AIC and AICc. Under H , there is a satisfactory weak selection of M with a rather important representation of M 1 base

Figure 1

 1 Figure 1 Decision path used to choose the best model for the selection strategies based on tests.

Figure 2 Figure 3

 23 Figure 2 Percentage of data sets simulated under H where the model with no gene effect (M ) is not selected versus the same percentage under H 1 base 0 for the eight model selection strategies using FOCE for N 40. The vertical line corresponds to a value of 5 and the horizontal line to a = %

  CT TT), complete model.mult {β 0 = β 1 ≠ β 2 ≠ } ≠ ≠

  CT TT), complete model.mult {β 0 = β 1 ≠ β 2 ≠ } ≠ ≠

  different model building strategies described above based either on tests or selection criteria are evaluated on data sets with the N 40 design. Results are reported as the percentage of data sets for which each model of the exon 26 polymorphism effect is selected, and

	=
	this with simulation under H or under H . The correct model is M 0 1	base	under H and M 0	mult	under H . The results of the strategies obtained 1
	by simulation under H are not used to modify the strategies under H as there are no simple correction methods for over-selecting 0 1
	covariate models under H .

0

Table 1

 1 Type I error for each test and for each algorithm.K is the number of data sets on which the test could be performed. Power for each test for N 40. K is the number of data sets on which the test could be performed. The corrected power was obtained using the fifth percentile of the empirical distribution of the test statistic under H as the cut-off value for the test.

		N 40 =
	Test	Algorithm

* Estimate significantly different from 5 . % Table 2 % Corrected power ( ) 0

Table 3

 3 Percentage of data sets simulated under H with a design of N 40 for which each model is selected.

	0	=
	Method	Algorithm	K

Table 4

 4 Percentage of data sets simulated under H with N 40 for which each model is selected. K is the number of data sets on which the test could be performed. 1 M :1 (CC CT TT) model with no gene effect.

	1	=
	Method	Algorithm	K

13

13

13

13

13

13

13

13