PhD Fokko J Bosker
email: f.j.bosker@psy.umcg.nl

PhD Catharina A Hartman

PhD Ilja M Nolte

Bram Prins

PhD Danielle Posthuma

PhD Tineke Van Veen

PhD Gonneke Willemsen

Roel De Rijk

PhD Eco J De Geus

Jg Witte

Hoogendijk

PhD Peter Terpstra

PhD Roel H Derijk

MD Witte J Hoogendijk

MD Patrick F Sullivan

PhD Brenda W Penninx

PhD Harold Snieder

MD Willem A Nolen

Poor replication of candidate genes for major depressive disorder using genome-wide association data

Keywords: candidate genes, genome-wide association study, major depressive disorder, replication, single nucleotide polymorphisms

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

Introduction

Major depressive disorder (MDD) is a multi-factorial disease, with both genetic and environmental factors likely to play a role in its etiology. The Netherlands Study of Depression and Anxiety (NESDA; www.nesda.nl) and the Netherlands Twin Registry (NTR; www.tweelingenregister.org) took part in the Genetic Association Information Network (GAIN; http://www.fnih.org/GAIN) project to enable a genome-wide association study (GWAS) using a 600K Perlegen chip [START_REF] Manolio | New models of collaboration in genome-wide association studies: the Genetic Association Information Network[END_REF] . Within the GAIN MDD study [START_REF] Boomsma | Genomewide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects[END_REF] 1862 participants with a diagnosis of MDD and 1860 controls at low liability for MDD were selected for genome-wide genotyping.

A GWAS approach allows a hypothesis-free search for potential new susceptibility genes. The downside of a GWAS is that a strict statistical adjustment for the large number of single nucleotide polymorphisms (SNPs; in the GAIN MDD study 435,291 SNPs) is required before associations can be considered significant on a genome-wide level [START_REF]The Wellcome Trust Case-Control Consortium: Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls[END_REF] and replication of such findings in independent cohorts is mandatory to exclude false positive findings [START_REF] Sullivan | Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo[END_REF] .

Another potential use of a GWAS is to use the results for a large-scale replication study of previous candidate gene studies. Application of such prior knowledge within the context of a GWAS allows less stringent significance thresholds than those for the hypothesis-free GWAS analysis [START_REF] Arking | A common genetic variant in the NOS1 regulator NOS1AP modulates cardiac repolarization[END_REF][START_REF] Roeder | Using linkage genome scans to improve power of association in genome scans[END_REF] . So far, the role of candidate genes in MDD has been the subject of many association studies. Unfortunately, there is little consistency between candidate gene studies for multifactorial diseases such as MDD (e.g., [START_REF] Hirschhorn | A comprehensive review of genetic association studies[END_REF][START_REF] Munafo | Candidate gene studies in the 21st century: meta-analysis, mediation, moderation[END_REF] ).

In the current study we attempted to replicate significant findings from previous candidate gene studies in MDD. To this end we conducted a systematic review of the literature and selected those genes that were reported to be significantly associated with MDD at least once. The GWAS data from the GAIN MDD study were used to screen all the identified candidate genes in two ways: 1) for association with the specific SNPs reported in the literature; and 2) for association with any of the SNPs (genotyped or imputed) from the Perlegen chip within the identified genes.

Materials and methods

Selection of studies reporting on candidate genes

In order to prevent any bias from the results of the GAIN MDD study, we identified candidate genes for MDD before the results from the GAIN MDD study became available. Therefore, a so called "enhanced search" was performed in Medline® via Pubmed on September 1, 2007 using the following search terms: (("genes" [ hits. We additionally scrutinized the reference list of the meta-analysis of genetic studies on MDD by Lopez-Leon and colleagues [START_REF] Lopez-Leon | Meta-analyses of genetic studies on major depressive disorder[END_REF] that appeared on-line on 16 October, 2007, shortly after the end date of our search, resulting in an additional 39 hits of possibly relevant papers.

These authors used somewhat broader search terms than we did, and their search ran until June 2007; Therefore, as a final check, we searched the literature using their search terms for the period June 2007-September 2007 not covered in their paper. This yielded an additional 110 hits. Of all these papers we retrieved the abstracts, and if considered relevant, the full paper.

In the next step, we selected all papers fulfilling the following five inclusion criteria: 1) the study had to be a candidate gene case-control association study. Linkage and finemapping studies were excluded; 2) the primary diagnosis of the patients in the candidate gene study had to be MDD to enable comparison with the GAIN MDD study. Therefore, we excluded studies involving patients (i) with a depressive episode in the course of bipolar disorder, (ii) with a primary psychotic disorder such as schizophrenia and a secondary depression, (iii) with a seasonal affective disorder not being MDD with a seasonal pattern, (iv) with a primary anxiety disorder (i.e. panic disorder, agoraphobia or social phobia) or obsessive compulsive disorder and a secondary depression, and (v) with MDD plus an additional specific comorbid disorder or condition, e.g. MDD plus alcoholism. However, we allowed subgroups within MDD, e.g. MDD in women or in men, recurrent MDD, or early onset MDD; 3) the sample of a specific study consisted of at least 30 patients with MDD and 30 healthy controls; 4) the findings on the association with MDD of any variant within the candidate gene (either a SNP, a microsatellite marker or a haplotype) had reached a statistical significance at the level of p < 0.05; and 5) the genetic association had to be with the diagnosis of MDD and not with other aspects such as associated personality features (e.g. neuroticism) or factors related to treatment response. This resulted in 78 papers.

In order not to miss potential true positive findings, we did not exclude candidate genes with associated markers that had low allele frequencies or that deviated from Hardy-Weinberg equilibrium. Two investigators (FB, CH) independently made a selection from the initial list of papers using the abovementioned criteria. When both reached consensus, the paper was included or excluded; in case of disagreement, consensus was obtained with assistance of two other authors (WN, HS). As a final step, one author (IN) double checked whether all selected papers fulfilled the aforementioned selection criteria 1 through 5. Figure 1 shows a flowchart of how we retrieved the 78 papers for the present study.

Bioinformatic tools

For many SNPs no reference SNP identification number (rs-id) was given in the original papers, but codes based on position (e.g., 677C/T in MTHFR or Tyr129Ser in HTR3B) or even own codes (e.g., s1-s5 in AVPR1B). In order to retrieve rs-ids for these SNPs, we used searches in PubMed or in NCBI's dbSNP database (http://www.ncbi.nlm.nih.gov/SNP/) using the "Geneview" option in conjunction with NCBI's MapViewer with the human genome assembly build 37.1 (http://www.ncbi.nlm.nih.gov/projects/mapview/map). Six SNPs remained that could not be easily found in this way. We contacted the corresponding authors of the papers and used NCBI's Primer-Blast and SNP-Blast (http://blast.ncbi.nlm.nih.gov/Blast.cgi) to map these SNPs using the provided primer sequences.

Sample

The 1,862 MDD cases included in the GAIN MDD study were mainly from NESDA, a longitudinal cohort study designed to be representative of individuals with depressive and/or anxiety disorders [START_REF] Penninx | The Netherlands Study of Depression and Anxiety (NESDA): Rationale, Objectives and Methods[END_REF] and were recruited from mental health care organizations (N=785), primary care (N=603), and community samples (N=314). Additional cases came from the NTR (N=160). Regardless of recruitment setting, similar inclusion and exclusion criteria were used to select MDD cases. Inclusion criteria were a lifetime diagnosis of MDD according to DSM-IV [START_REF]American Psychiatric Association: Diagnostic and Statistical Manual of Mental Disorders[END_REF] as diagnosed via the Composite International Diagnostic Interview (CIDI Version 2.1. [START_REF]World Health Organization: Composite International Diagnostic Interview (CIDI), Version 2.1[END_REF] ), age 18-65 years, and self-reported western European ancestry. Persons who were not fluent in Dutch and those with a primary diagnosis of a psychotic disorder, obsessive compulsive disorder, bipolar disorder, or severe alcohol or substance use disorder were excluded.

Most of the 1,860 control subjects were from the NTR (N=1,703) and additional controls from NESDA (N=157). Longitudinal phenotyping in NTR includes assessment of depressive symptoms (via multiple instruments), anxiety and neuroticism. Inclusion required no report of MDD at any measurement occasion and never scoring high (>0.5 SD) on a general factor score based on combined measures of neuroticism, anxiety and depressive symptoms. Finally, controls and their parents were required to have been born in the Netherlands or Western Europe. Only one control per family was selected. NESDA controls had no lifetime diagnosis of MDD or an anxiety disorder as assessed by the CIDI and reported low depressive symptoms at baseline. For more details see [START_REF] Boomsma | Genomewide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects[END_REF] .

Genotyping

Perlegen Sciences (Mountain View, CA, USA) performed all genotyping according to strict standard operating procedures. DNA samples from cases and controls were randomly assigned to plates, shipped to Perlegen and identified only by barcode. High-density oligonucleotide arrays were used yielding 599,164 SNPs. Eight SNPs with duplicate numbers were deleted and 73 mitochondrial SNPs were removed for later analysis. From the remaining 599,083 SNPs on the Perlegen chip 435,291 passed quality control. Two hundred eighty subjects were excluded because of various quality control issues leaving 1,738 cases (93.4%) and 1,802 controls (96.9%) in the final analysis dataset. For more details see [START_REF] Boomsma | Genomewide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects[END_REF][START_REF] Sullivan | Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo[END_REF] .

Imputation

Not all SNPs selected from the literature were present on the genotyping array. Based on the linkage disequilibrium (LD) structure between SNPs we followed an imputation procedure to predict non-genotyped SNPs using the HapMap CEU data (release 22, build 36) as the reference database. The imputation was performed by IMPUTE version 0.3.2 using the default settings and the recommended number 11418 for the effective population size of Caucasians [START_REF] Marchini | A new multipoint method for genomewide association studies by imputation of genotypes[END_REF] . In this way we extended the genome-wide autosomal SNP dataset from 427,049 to 2,467,430 SNPs. For our candidate genes this meant an extension from 851 to 4,955 SNPs. However, the quality of the imputation was low for 85 SNPs (SNPTEST proper_info <0.5). These SNPs were excluded leaving 4,870 SNPs for analysis. No SNP had a minor allele frequency <1%.

Association test

Association between MDD and the autosomal SNP data was tested by means of a frequentist case-control test provided in the software package SNPTEST version 1.1.4 using allele dosages with sex as a covariate in order to adjust for the slight imbalance in percentage of females between cases (69.6%) and controls (62%) [START_REF] Boomsma | Genomewide association of major depression: description of samples for the GAIN Major Depressive Disorder Study: NTR and NESDA biobank projects[END_REF][START_REF] Sullivan | Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo[END_REF] and the "proper" option to account for the uncertainty of the genotypes that were imputed [START_REF] Marchini | A new multipoint method for genomewide association studies by imputation of genotypes[END_REF] . In addition, 7,988 genotyped SNPs on the X chromosome were analyzed in PLINK version 1.03 [START_REF] Purcell | PLINK: a tool set for whole-genome association and population-based linkage analyses[END_REF] by means of logistic regression with sex as a covariate. SNPs on the Y chromosome (n=15) and SNPs mapped to ambiguous locations (n=239) were not analyzed.

Odds ratios (ORs) and 95% confidence intervals (CIs) were calculated for the significant candidate SNPs according to the disease model from the original article using counted or estimated numbers for genotyped or imputed SNPs, respectively. In order to establish true replication, we checked whether the effect was for the same allele and in the same direction.

Permutation procedure

In order to facilitate interpretation of the significance of our findings for each SNP we calculated three p-values by permutation: (i) a crude uncorrected significance; (ii) a genewide significance (corrected for all SNPs in the gene), and (iii) an overall significance (corrected for all SNPs in all selected genes). The crude uncorrected p-value was determined as the fraction of permutations that yielded a p-value that was smaller than the observed one. This p-value was determined in order to validate the permutation procedure since it should be similar to the p-value calculated by SNPTEST or PLINK. The gene-wide significance of a particular SNP was defined as the fraction of permutations that the p-value of any SNP in the candidate gene concerned was smaller than the one observed for that SNP. This procedure corrects the p-value for testing multiple SNPs within a gene. Based on the rationale that each selected gene is a candidate for MDD (hypothesis-driven) no correction for all SNPs in the entire gene set is necessary. Nonetheless, we also calculated an overall significance (i.e., corrected for all SNPs in all candidate genes) by determining an overall p-value for each SNP as the fraction of permutations for which any of the SNPs in any of the candidate genes had a p-value smaller than the observed one. For each of the three pvalues we conducted 10,000 permutations. Case-control status was randomly permuted within males and females separately, hence leaving the number of affected males and females intact.

In addition to the three SNP specific p-values, we also determined whether the number of SNPs within a candidate gene with an original p-value <0.05 based on SNPTEST or PLINK (see above) was higher than expected. For each permutation the number of SNPs NESDA-NTR GAIN Candidate gene replication study 9 in each candidate gene that was significant at 0.05 was recorded. The fraction of permutations with a higher number of significant SNPs than originally observed determined the significance of the number of significant SNPs of that candidate gene. By permuting case-control status the LD structure of the genes is preserved, that is, the resulting significance is corrected for possible high correlation between SNPs. The number of permutations for this test was 10,000 as well.

Results

Literature search

The 78 papers which resulted from our systematic literature search reported 115 statistically significant differences between MDD cases and healthy controls in 57 genes: for 96 SNPs, 7 microsatellite markers (i.e., length polymorphisms), 11 haplotypes, and one protein polymorphism in the haptoglobin gene [START_REF] Maes | Haptoglobin phenotypes and gene frequencies in unipolar major depression[END_REF] (Table 1). Twenty-nine SNPs were identified by an rs-id, while 67 SNPs were only specified by a location code, restriction enzyme or author designed code. For 64 of these 67 SNPs, rs-ids could be obtained, while two rs-ids were not found (246G/A in CCKAR [START_REF] Koks | Analysis of SNP profiles in patients with major depressive disorder[END_REF] and -7054C/A in DRD2 [START_REF] Koks | Analysis of SNP profiles in patients with major depressive disorder[END_REF] ). In addition, for SNP 1463G/A in TPH2, others have tried to replicate the association, but the SNP appeared to be nonexistent [START_REF] Zhang | Loss-offunction mutation in tryptophan hydroxylase-2 identified in unipolar major depression[END_REF][START_REF] Blakely | Overview: a rare opportunity or just one less reason to be depressed[END_REF] . Thus, we could map 93 SNPs.

For the 7 microsatellite markers no LD data with SNPs are known and hence we cannot determine whether these are covered by any of the available SNPs in the corresponding genes. Among them is the 5-HTTLPR 44 base pair deletion in the promoter region of the serotonin transporter gene SLC6A4 (SERT), which has attracted considerable attention in various previous candidate gene studies for MDD.

Association in the GAIN MDD data with specific SNPs from the literature

Sixty one of the 93 selected SNPs in the 57 candidate genes were either present (n= 18) on the Perlegen array or could be imputed (n= 43). Four additional SNPs could be tagged by other available SNPs (Table 2). Two of these were not genotyped in the CEU population of the HapMap Phase2 (release 22, build 36) project, but were available in the JPT and CHB populations in which they showed high LD (r 2 =0.97 and 1.00, respectively) with at least one other available SNP. The two other SNPs were genotyped in the CEU population and could be tagged by available SNPs with r 2 =1, but were for unreported reasons not included in the HapMap reference file used for the imputation procedure as provided on the website of the imputation software package IMPUTE (http://www.stats.ox.ac.uk/~marchini/software/gwas/impute.html). Twenty eight SNPs were neither genotyped directly nor imputed nor could be tagged, hence were not available for replication.

We investigated which of the 65 available or tagged SNPs could be confirmed in our data (Table 2). SNPs rs12520799 in C5orf20 (DCNP1) (p=0.038; OR AT=1.10, 95% CI 0.95-1.29; OR TT=1.21, 95% CI 1.01-1.47) and rs16139 in NPY (p=0.034; OR C allele=0.73, 95% CI 0.55-0.97) replicated the reported effect in the literature. In addition, three SNPs in the ACE gene selected from the literature, which were in strong LD with each other (r 2 =0.78-1.00), were also significant in the GAIN MDD sample (rs4333: p=0.029; rs4329: p=0.030; rs4461142: p=0.036). However these three SNPs showed effects in the opposite direction compared to previous results [START_REF] Baghai | Polymorphisms in the angiotensin-converting enzyme gene are associated with unipolar depression, ACE activity and hypercortisolism[END_REF] (Baghai, personal communication). Hence these SNPs do not represent a true replication.

Candidate genes from the literature in the GAIN MDD study

As in different populations different SNPs might play a role, we also studied all SNPs present on the genotyping array or available through imputation in an area from 5kb upstream to 5kb downstream of each selected candidate gene in order to cover the promoter and 3' untranslated region, respectively. Forty nine of the 57 genes were covered by one or more SNPs that were present on the Perlegen chip. For another six genes no genotyped SNPs were available but imputed ones were. Neither genotyped nor imputed SNPs were available for two genes (AVPR1B and CHRFAM7A). In total 4,870 SNPs, of which 820 were genotyped and 4,019 were imputed, covered the 55 candidate genes ranging from 1 SNP for AR to 642 SNPs for PDE11A (Figure 2). We noted a significant correlation of 0.75 between gene size and the number of genotyped SNPs (including genes on the X chromosome) and a correlation of 0.96 between gene size and total number of SNPs (excluding genes on the X chromosome, since SNPs on the X chromosome were not imputed).

For 28 of the 55 genes one or more SNPs were found to be different between MDD cases and healthy controls in the GAIN MDD GWA scan at a significance level of p<0.05 (Table 3). The remaining 27 genes were not associated with MDD as none of the SNPs reached p<0.05. The smallest p-value was observed for SNP rs769178 in the TNF gene region (p=0.00029; OR T allele =1.35, 95%CI 1.13-1.63). The minor allele T was observed in 8.0% of the MDD cases and 6.1% of the controls.

With genes covered by a large number of SNP markers, the expected minimal pvalue will decrease purely as a result of chance alone as a function of the number of SNPs. Thus, we used gene-wide p-values from the permutation procedure that corrects for the number of SNPs within a candidate gene and only SNP rs769178 in TNF remained significant (p=0.0034) (Figure 3; Table 3. The second strongest associated SNP, which was observed in DISC1 (rs7533169: p-value=0.0025), became non-significant after this correction (gene-wide p-value=0.28), because in this gene there were 491 SNPs and apparently the small p-value was observed purely based on chance. In addition to the gene-wide significance, we also determined overall significance per SNP accounting for testing 4,870 SNPs in 55 candidate genes. In that case none of the SNPs remained significant, not even the TNF SNP rs769178 (overall p=0.33).

The significance levels per SNP are one way of testing the true value of the selected candidate genes. If many SNPs in a candidate gene are associated, this could also indicate potential involvement of a gene in the disorder under study [START_REF] Arking | Multiple independent genetic factors at NOS1AP modulate the QT interval in a multi-ethnic population[END_REF] . We noticed for instance that in the norepinephrine transporter SLC6A2 (NET) 19 of the 100 SNPs were significant at a level of 0.05 (Table 3). The permutation procedure that tested whether the number of SNPs that were significant at 0.05 was larger than expected, revealed that this finding for SLC6A2 was NESDA-NTR GAIN Candidate gene replication study 12 indeed significant (p=0.039). Figure 4 shows that the 19 SNPs are scattered in the right half of the gene. Seventeen of them lie in three independent haplotype blocks and the other 2 SNPs are not in LD with the haplotype blocks or with each other. For none of the other candidate genes this permutation test revealed a significant result.

Discussion

Several genome-wide linkage studies of MDD have been published (reviewed in 2 ) but the GAIN MDD study is among the first genome-wide association studies (GWAS) in MDD [START_REF] Sullivan | Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo[END_REF][START_REF] Muglia | Genome-wide association study of recurrent major depressive disorder in two European case-control cohorts[END_REF] .

We used these GWAS data as a large scale replication of previously reported significant findings from candidate gene studies in MDD. To this end, we first conducted a systematic review of the literature and identified a total of 57 genes for which a significant association with MDD has been reported at least once. Fifty-five of these genes could be included in our replication study, with either genotyped or imputed SNPs available from the GWAS data.

With a sample size (1,738 cases and 1,802 controls ) by far exceeding all previous candidate gene studies (with a mean sample size of 164 MDD cases and 252 controls, and only 1 study with > 1000 MDD cases and controls) [START_REF] Lucae | P2RX7, a gene coding for a purinergic ligand-gated ion channel, is associated with major depressive disorder[END_REF] , we found minimal support for involvement in MDD for only 4 out of 55 candidate genes: C5orf20 (p=0.038), NPY (p=0.034), TNF (p=0.0034), and SLC6A2 (p=0.039). Replication of these genes was based on three different statistical approaches. First, the involvement of C5orf20 (rs12520799) and NPY (rs16139) constituted a direct replication at SNP level of previously identified SNPs associated with MDD. Second, studying the selected candidate genes for all SNPs present on the genotyping array, TNF (rs76917) was identified as the only gene for which the association with MDD remained gene-wide significant. Third, the potential involvement of SLC6A2 (NET) in MDD was derived on the basis of a statistically significant number of SNPs (i.e., 19 out of 100) associated with MDD within this gene, which could not be explained solely by high LD between the SNPs. Note that the prior evidence in the selected studies for the involvement of C5orf20, NPY, TNF, and SLC6A2 in MDD did not stand out as particularly strong, i.e. odds ratios were mostly in the moderate range (minimum OR = 0.33; maximum OR = 2.41), and none of the p-values in these studies was smaller than p<0.001, not even for the largest study on C5orf20 with a total sample size of 864. NPY was the one exception where the allele was present in 6.3 % of the controls but not in the patients, thus suggesting a strong effect (OR=0). However, this was based on only 51 patients [START_REF] Willis-Owen | DCNP1: a novel candidate gene for major depression[END_REF][START_REF] Heilig | Decreased cerebrospinal fluid neuropeptide Y (NPY) in patients with treatment refractory unipolar major depression: preliminary evidence for association with preproNPY gene polymorphism[END_REF][START_REF] Inoue | Positive association between T-182C polymorphism in the norepinephrine transporter gene and susceptibility to major depressive disorder in a japanese population[END_REF][START_REF] Jun | Possible association between -G308A tumour necrosis factor-alpha gene polymorphism and major depressive disorder in the Korean population[END_REF] .

An important question is why so few candidate genes were replicated by our GWAS, while the sample size of our study was so much larger than any of the 78 selected candidate gene studies. One possible -and in our opinion most likely -explanation is publication bias: positive results have a better chance of being published than negative results. This would also imply that many previously reported positive findings were actually false positive findings (type 1 errors), probably amplified by insufficient correction for multiple testing. Testing this, e.g. with funnel plots, is not appropriate for the approach followed here as our literature search strategy did not include negative candidate gene studies. Second, given that MDD is a rather heterogeneous disorder, and that it was diagnosed with different instruments across studies, previous and current study samples may have differed phenotypically (see [START_REF] Hek | The PCLO gene and depressive disorders: replication in a population-based study[END_REF] for a discussion of this point). A third explanation is that associations between genes and etiologically complex diseases may depend on genetic (gene*gene interactions or epistasis) and environmental (gene*environment interactions) contexts, which may differ in samples from different populations [START_REF] Martinez | Gene-environment interactions in asthma: with apologies to William of Ockham[END_REF] . Thus, previous and current samples may have been different genetically or in their contextual factors.

Given the modest support for the replicated candidate genes (one p-value <0.01, three p-values <0.05) it is possible that even the four replications of our analysis are false positive results. With the number of SNPs (n=65) been tested and under the null hypothesis of no true associations of any of the candidate SNPs in previous studies, one would expect three significant findings. Thus, it is possible that our two significant findings in C5orf20 and NPY are false positive results. In addition, under a similar null hypothesis for the 55 candidate genes, two significant findings were expected implying that the single gene-wide significant result of TNF might also be a false positive result. This is supported by the lack of overall significance for the TNF SNP (p=0.33). And finally, also the significant finding for NESDA-NTR GAIN Candidate gene replication study 14 SLC6A2, the only gene showing a larger number of significant SNPs than expected, might be a chance finding.

On the other hand, the above calculations may be too conservative when assuming that at least some of the previously found associations of candidate genes with MDD were true findings. In that case, our approach of replicating candidate genes is more or less hypothesis-driven, thus not requiring the same multiple testing penalty as the genome-wide approach [START_REF] Sullivan | Genome-wide association for major depressive disorder: a possible role for the presynaptic protein piccolo[END_REF] . Nonetheless, the few replications in our study are rather sobering and in order to uncover whether the replicated SNPs and genes are truly associated to MDD confirmation in independent samples is crucial. As such, meta-analytical results from the Psychiatric GWAS Consortium are also eagerly awaited [START_REF]A framework for interpreting genomewide association studies of psychiatric disorders[END_REF] .

From the most recent meta-analysis of genetic studies on MDD, Lopez-Leon et al. [START_REF] Lopez-Leon | Meta-analyses of genetic studies on major depressive disorder[END_REF] concluded that statistically significant evidence exists for six MDD susceptibility genes, i.e., APOE, DRD4, GNB3, MTHFR, SLC6A3, and SLC6A4. Our study offers little support for these genes. Given our sample size we had 80% power to detect odds ratios of >1.15 for allele frequencies >5% under an additive disease model. All of the above genes meet these criteria. The association of GNB3 and MTHFR with MDD was directly tested but could not be replicated in the GAIN MDD sample, although the effect of MTHFR showed a trend in the expected direction (p=0.074; OR TT versus CC=1.14). In addition, APOE, DRD4, and SLC6A4 were not associated with MDD in our study, but the previously identified genetic variants were length polymorphisms instead of SNPs. Hence these could not directly be tested and since it is unsure whether these variants are tagged by the SNPs in the corresponding genes we cannot refute the associations. SLC6A3 was not identified in our literature search as a candidate gene for MDD, because individual studies did not report significant effects for this gene [START_REF] Frisch | Association of unipolar major depressive disorder with genes of the serotonergic and dopaminergic pathways[END_REF][START_REF] Manki | Dopamine D2, D3 and D4 receptor and transporter gene polymorphisms and mood disorders[END_REF] and it only reached significance in the pooled metaanalysis. Post-hoc we also tested for significance in our sample and we did not find any association (59 SNPs, most significant SNP: rs27072, P=0.096).

In the context of the non-replication of many of the selected candidate genes as discussed earlier, two limitations of our study need to be noted. First, we did not have direct NESDA-NTR GAIN Candidate gene replication study 15 or indirect information on one third of the candidate SNPs as these were not present on the genotyping chip and could not be imputed using the HapMap CEU data. Second, we were unable to test candidate length polymorphisms previously associated with MDD. Among these was the 44 base pair insertion/deletion polymorphism (or 5-HTTLPR) in the promoter region of the serotonin transporter gene SLC6A4 (SERT). Length polymorphisms are often difficult to tag with single SNPs because LD information between SNPs and length polymorphisms is either unavailable or LD with SNPs is insufficiently strong [START_REF] Su | Serotonin Transporter Gene, Depressive Symptoms and Interleukin-6[END_REF] . However, a recent study by Wray and colleagues [START_REF] Wray | Accurate, Large-Scale Genotyping of 5HTTLPR and Flanking Single Nucleotide Polymorphisms in an Association Study of Depression, Anxiety, and Personality Measures[END_REF] identified a two-SNP haplotype proxy for 5HTTLPR.

In conclusion, the GWAS data of the GAIN MDD study failed to replicate all but four of the previously reported candidate gene associations with MDD. However, given the number of candidate SNPs and genes that were tested, even these significant may well be false positives, implying that we found no replication at all. This might point to publication bias and false-positive findings in previous studies, but also to heterogeneity of the MDD phenotype as well as variations in contextual genetic or environmental factors.

  TIAB] NOT Medline[SB]) OR "genes"[MeSH Terms] OR gene[Text Word]) OR ("genes"[MeSH Terms] OR genes[Text Word]) OR snp[All Fields] OR ("single nucleotide polymorphism"[Text Word] OR "polymorphism, single nucleotide"[MeSH Terms] OR snps[Text Word]) OR ("genetic polymorphism"[Text Word] OR "polymorphism, genetic"[MeSH Terms] OR polymorphism[Text Word]) OR polymorphisms[All Fields] AND ("major depressive disorder"[Text Word] OR "depressive disorder, major"[MeSH Terms] OR major depression[Text Word]) AND "humans"[MeSH Terms]. This resulted in 641

Figure 1 :

 1 Figure 1: Selection procedure of studies reporting on candidate genes. PM, PUBMED; LL, Lopez-Leon.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Coverage by genotyped (black bar) and imputed (grey bar) SNPs and size (white circles) of the 57 selected candidate genes (above: genes with 50 or less SNPs, below: genes with more than 50 SNPs) Figure 3: Observed minimal p-value versus the expected minimal p-value based on the number of SNPs in the gene

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Coverage by genotyped (black bar) and imputed (grey bar) SNPs and size (white circles) of the 57 selected candidate genes (above: genes with 50 or less SNPs, below: genes with more than 50 SNPs)

Table 1 : Candidate genes and associated polymorphisms or haplotypes based on the systematic literature search. Sample sizes (N cases/N controls) as well as p-values for allelic (P allele), genotypic (P genotype) and/or multimarker tests (P multimarker) are shown. Odds ratios (ORs) are given whenever possible.

 1 SNP is located more than 5kb away from gene and is not included in our analyses b OR was not reported in corresponding article, but could be estimated from genotype counts or frequencies c allele was not observed in cases, which resulted in OR = 0.00

	Gene	HUGO name	Chromosome position	Gene size (kb)	Variant	rs-id	Refs	N cases	N controls	P allele	OR allele	P genot.	OR genot.		P multi-marker	OR
	MTHFR	MTHFR	1p36.3	20.3	677C/T	rs1801133	34	32	419	0.005	2.8
							35	100	89				0.03	1.90
	5-HT1A	HTR1A	5q11.2-q13	1.3	-1019C/G	rs6295	36	129	134	0.0006	1.84	b	0.0017	het:1.54 hom:3.75	b
	5-HT2A	HTR2A	13q14-q21	62.7	-1438A/G	rs6311	37	71	157	0.018	0.62	0.039	0.47
							38	189	148	0.007	1.52	0.005
					102T/C	rs6313	39	120	131	0.03			0.01	2.4
	5-HT2C	HTR2C	Xq24	326.1 68G/C = Cys23Ser	rs6318	16	177	160	0.02	1.84 b
							40	513	901	0.006	n.a.
	5-HT3B	HTR3B	11q23.1	41.7	Tyr129Ser (3B06)	rs1176744	41	47	90				0.027	n.a.
					Ala154Ala (3B08)	rs2276305	41	47	90	0.015	2.27	0.011	n.a.
					IVS6+72A>G (3B10)	rs2276307	41	47	90				0.023	n.a.
					IVS6+161C>G (3B11)	rs2276308	41	47	90				0.031	n.a.
					3B04+3B06+3B08+3B09 +3B10+3B11+3B12	rs1176746+rs1176744 +rs2276305+rs3782025	41	47	90							0.002	2.11
	ACE	ACE	17q23.3	20.5	dbSNP4295	rs4295	19	642	608				0.0011	rec:1.87	b
							42	92	92	<0.0001			0.001			0.0001 n.a.
					dbSNP4309	rs4309	42	92	92	0.03	1.59

a

Table 3 : Gene-wide results for the 57 candidate genes identified by the systematic literature search

 3 -values were determined using SNPTest for the autosomes and logistic regression in PLINK for the X chromosome. Permutations (N=10,000) were performed to determine uncorrected SNP p-values, gene-wide corrected p-values, overall corrected p-values and p-values for number of SNPs with p<0.05. All p-values were calculated assuming an additive model. P-values <0.05 are shown in bold.

					minimal p	minimal p				#SNPs at	p-value
	gene	#genotyped #imputed	#total	(analytical)	(permuted) gene-wide p	overall p	rs-id	p<0.05	#SNPs
	TNF	2	18	20	0.00029	0.00094	0.0034	0.34	rs769178	1	0.32
	DISC1	105	386	491	0.0025	0.0015	0.29	0.95	rs7533169	18	0.56
	SLC6A2	25	75	100	0.0061	0.0058	0.19	1.00	rs5564	19	0.039
	DRD2	12	68	80	0.0079	0.0077	0.15	1.00	rs7131056	3	0.32
	HTR3B	6	33	39	0.0081	0.0077	0.11	1.00	rs7945926	3	0.21
	CHRM2	43	157	200	0.0097	0.0081	0.32	1.00	rs6944132	3	0.54
	MTHFR	6	37	43	0.0098	0.0080	0.11	1.00	rs3737967	4	0.18
	PDE10A	107	392	499	0.013	0.014	0.74	1.00	rs7741623	8	0.80
	APOE	0	5	5	0.015	0.012	0.057	1.00	rs10119	1	0.11
	TACR1	31	141	172	0.016	0.019	0.50	1.00	rs1012686	24	0.10
	TPH2	21	134	155	0.017	0.016	0.29	1.00	rs17722134	6	0.29
	PDE11A	66	576	642	0.018	0.016	0.63	1.00	rs6433706	11	0.57
	DRD3	13	41	54	0.018	0.021	0.23	1.00	rs1800827	3	0.27
	PDE5A	7	108	115	0.019	0.020	0.20	1.00	rs13104219	21	0.08
	NR3C1	14	60	74	0.024	0.025	0.37	1.00	rs4607376	12	0.11
	MAOA	5	0	5	0.026	0.028	0.078	1.00	rs5906883	1	0.11
	ESR1	50	272	322	0.027	0.025	0.84	1.00	rs17081749	1	0.93
	ACE	3	28	31	0.028	0.031	0.21	1.00	rs4344	8	0.08
	NPY	6	38	44	0.030	0.031	0.29	1.00	rs1073317	6	0.15
	PDE9A	25	133	158	0.033	0.031	0.69	1.00	rs13046735	1	0.75
	HTR2A	35	100	135	0.037	0.037	0.72	1.00	rs2770292	3	0.56
	C5orf20	10	35	45	0.038	0.038	0.36	1.00	rs12520799	5	0.13
	TPH1	6	24	30	0.040	0.038	0.25	1.00	rs10741734	2	0.19
	PLA2G4A	26	92	118	0.040	0.042	0.68	1.00	rs7555140	1	0.68
	BCR	14	110	124	0.041	0.045	0.69	1.00	rs11090231	1	0.67
	ESR2	43	155	198	0.042	0.040	0.80	1.00	rs1255949	2	0.68
	P2RX7	21	87	108	0.047	0.049	0.62	1.00	rs11065464	1	0.57
	AGTR1	8	44	52	0.054	0.049	0.58	1.00	rs12721331	0	0.77
	PDE6C	18	52	70	0.055	0.056	0.71	1.00	rs10882288	0	0.84
	CYP2C9	8	47	55	0.057	0.056	0.48	1.00	rs1799853	0	0.72
	CCKAR	4	17	21	0.060	0.058	0.33	1.00	rs11726159	0	0.64

  [START_REF] Penninx | The Netherlands Study of Depression and Anxiety (NESDA): Rationale, Objectives and Methods[END_REF] (p-value) on the y-axis. Nineteen SNPs have a p-value <0.05, which is significantly more than expected (p=0.039). Seventeen of them lie in three haplotype blocks (squares, upward triangles [light grey], and downward triangles [dark grey]) (r 2 values HapMap CEU: squares-upward triangles: 0.020, squaresdownward triangles: 0.661, and upward triangles-downward triangles: 0.013). The other 2 (black diamonds) are not in LD with the haplotype blocks or with each other. The light grey line shows the recombination rate in this area (axis on the right).

	641 articles identified from	149 articles identified through search
	PUBMED systematic search	terms and reference list of Lopez-Leon 9
	790 articles in total
	Reached Consensus
	FB/CH
	Yes	No
		41 articles
		(27 PM/14 LL)
		Additional Consensus
		WN/HS
		Yes
	74 articles	10 articles
	(67 PM/7 LL)	(8 PM/2 LL)
	84 articles in total
	Final check	Excluded 6
	IN		articles
	Included 78 articles

Figure 1

Selection procedure of studies reporting on candidate genes. PM, PUBMED; LL, Lopez-Leon.

Table 1 : Candidate genes and associated polymorphisms or haplotypes based on the systematic literature search. Sample sizes (N cases/N controls) as well as p-values for allelic (P allele), genotypic (P genotype) and/or multimarker tests (P multimarker) are shown. Odds ratios (ORs) are given whenever possible.

 1

	Gene	HUGO name	Chromosome position	Gene size (kb)	Variant	rs-id	Refs	N cases	N controls	P allele	OR allele	P genot.	OR genot.		P multi-marker	OR
	MTHFR	MTHFR	1p36.3	20.3	677C/T	rs1801133	34	32	419	0.005	2.8
							35	100	89				0.03	1.90
	5-HT1A	HTR1A	5q11.2-q13	1.3	-1019C/G	rs6295	36	129	134	0.0006	1.84	b	0.0017	het:1.54 hom:3.75	b
	5-HT2A	HTR2A	13q14-q21	62.7	-1438A/G	rs6311	37	71	157	0.018	0.62	0.039	0.47
							38	189	148	0.007	1.52	0.005
					102T/C	rs6313	39	120	131	0.03			0.01	2.4
	5-HT2C	HTR2C	Xq24	326.1 68G/C = Cys23Ser	rs6318	16	177	160	0.02	1.84 b
							40	513	901	0.006	n.a.
	5-HT3B	HTR3B	11q23.1	41.7	Tyr129Ser (3B06)	rs1176744	41	47	90				0.027	n.a.
					Ala154Ala (3B08)	rs2276305	41	47	90	0.015	2.27	0.011	n.a.
					IVS6+72A>G (3B10)	rs2276307	41	47	90				0.023	n.a.
					IVS6+161C>G (3B11)	rs2276308	41	47	90				0.031	n.a.
					3B04+3B06+3B08+3B09 +3B10+3B11+3B12	rs1176746+rs1176744 +rs2276305+rs3782025	41	47	90							0.002	2.11
	ACE	ACE	17q23.3	20.5	dbSNP4295	rs4295	19	642	608				0.0011	rec:1.87	b
							42	92	92	<0.0001			0.001			0.0001 n.a.
					dbSNP4309	rs4309	42	92	92	0.03	1.59

Table 3 : Gene-wide results for the 57 candidate genes identified by the systematic literature search

 3 -values were determined using SNPTest for the autosomes and logistic regression in PLINK for the X chromosome. Permutations (N=10,000) were performed to determine uncorrected SNP p-values, gene-wide corrected p-values, overall corrected p-values and p-values for number of SNPs with p<0.05. All p-values were calculated assuming an additive model. P-values <0.05 are shown in bold.

					minimal p	minimal p				#SNPs at	p-value
	gene	#genotyped #imputed	#total	(analytical)	(permuted) gene-wide p	overall p	rs-id	p<0.05	#SNPs
	TNF	2	18	20	0.00029	0.00094	0.0034	0.34	rs769178	1	0.32
	DISC1	105	386	491	0.0025	0.0015	0.29	0.95	rs7533169	18	0.56
	SLC6A2	25	75	100	0.0061	0.0058	0.19	1.00	rs5564	19	0.039
	DRD2	12	68	80	0.0079	0.0077	0.15	1.00	rs7131056	3	0.32
	HTR3B	6	33	39	0.0081	0.0077	0.11	1.00	rs7945926	3	0.21
	CHRM2	43	157	200	0.0097	0.0081	0.32	1.00	rs6944132	3	0.54
	MTHFR	6	37	43	0.0098	0.0080	0.11	1.00	rs3737967	4	0.18
	PDE10A	107	392	499	0.013	0.014	0.74	1.00	rs7741623	8	0.80
	APOE	0	5	5	0.015	0.012	0.057	1.00	rs10119	1	0.11
	TACR1	31	141	172	0.016	0.019	0.50	1.00	rs1012686	24	0.10
	TPH2	21	134	155	0.017	0.016	0.29	1.00	rs17722134	6	0.29
	PDE11A	66	576	642	0.018	0.016	0.63	1.00	rs6433706	11	0.57
	DRD3	13	41	54	0.018	0.021	0.23	1.00	rs1800827	3	0.27
	PDE5A	7	108	115	0.019	0.020	0.20	1.00	rs13104219	21	0.08
	NR3C1	14	60	74	0.024	0.025	0.37	1.00	rs4607376	12	0.11
	MAOA	5	0	5	0.026	0.028	0.078	1.00	rs5906883	1	0.11
	ESR1	50	272	322	0.027	0.025	0.84	1.00	rs17081749	1	0.93
	ACE	3	28	31	0.028	0.031	0.21	1.00	rs4344	8	0.08
	NPY	6	38	44	0.030	0.031	0.29	1.00	rs1073317	6	0.15
	PDE9A	25	133	158	0.033	0.031	0.69	1.00	rs13046735	1	0.75
	HTR2A	35	100	135	0.037	0.037	0.72	1.00	rs2770292	3	0.56
	C5orf20	10	35	45	0.038	0.038	0.36	1.00	rs12520799	5	0.13
	TPH1	6	24	30	0.040	0.038	0.25	1.00	rs10741734	2	0.19
	PLA2G4A	26	92	118	0.040	0.042	0.68	1.00	rs7555140	1	0.68
	BCR	14	110	124	0.041	0.045	0.69	1.00	rs11090231	1	0.67
	ESR2	43	155	198	0.042	0.040	0.80	1.00	rs1255949	2	0.68
	P2RX7	21	87	108	0.047	0.049	0.62	1.00	rs11065464	1	0.57
	AGTR1	8	44	52	0.054	0.049	0.58	1.00	rs12721331	0	0.77
	PDE6C	18	52	70	0.055	0.056	0.71	1.00	rs10882288	0	0.84
	CYP2C9	8	47	55	0.057	0.056	0.48	1.00	rs1799853	0	0.72
	CCKAR	4	17	21	0.060	0.058	0.33	1.00	rs11726159	0	0.64

tagged by rs12270070 (r 2 =1)

tagged by rs12270070 (r 2 =1

Acknowledgements:

The authors would like to thank Ting Wu for making the flow diagram in figure 1.

Conflict of interest: None


N Touati Moungla

L Létocart

A Nagih

An improving dynamic programming algorithm to solve the shortest path problem with time windows

Keywords:

An improving dynamic programming algorithm to solve the shortest path problem with time windows

Nora Touati Moungla, Lucas Létocart, Anass

Introduction

Dynamic Programming (DP) (1) is a well-established method to solve shortest path problems. In the classical shortest path problem, the Ford-Bellman algorithm assigns a single label to each node representing the cost which is gradually improved. An extension of this approach was proposed in [START_REF] Desrochers | An algorithm for the shortest path problem with resource constraints[END_REF] in the shortest path problem with ressource constraints context. For the shortest path problem with time windows, the routes must be compared in terms of arrival times at the nodes as well as costs: at each node i a label (arrival time at node i, cost of the route) is defined. Several labels have to be stored at node i to calculate the labels of other nodes according to the following optimality principle: If X sj is a minimum cost route from among all the routes from s to j arriving at j at time T j or before, and if (i, j) is its terminal arc, then the sub-route X si is a minimum cost route from among all routes from s to i arriving at i at time T jt ij or before, where t ij is the travel time of arc (i, j).

This procedure can leads to the computation of a large number of labels at each node of the graph. Dominance rules are then used to compare subroutes arriving at a same location for discarding some of them. Dynamic programming algorithms (2; 3) rely on labels extension and a dominance procedure [START_REF] Desrochers | A generalized permanent labelling algorithm for the shortest path problem with time windows[END_REF]. In label setting algorithms, nodes are treated once and labels associated to a treated node are kept until the end of the resolution process. In label correcting approaches, nodes are repeatedly treated and their labels extended to all feasible directions.

The dominance procedure applied on each new label requires in the worst case, the exploration of all efficient labels of the current node. Since the number of labels increases exponentialy as a function of the problem's size, dominance computing time increases too. We propose in this work a new variant of the dynamic programming algorithm for accelerating the dominance procedure by reducing the number of labels saved and by identifying efficiently dominated labels.

We present in section 2 the Shortest Path Problem with Time Windows (SPPTW) and describe in section 3 its resolution by the dynamic programming algorithm. We present in section 4 a new variant of this method called dynamic programming with blocs. We validate in section 5 this new approach by experiments on shortest path problem with time windows instances. In section 6 we conclude.

Problem description

Consider a connected graph G = (N , A), where N is the set of nodes and A is the set of arcs. With each arc (i, j) ∈ A is associated a real number cost c ij and a non negative integer duration t ij . Each node i ∈ N is characterized by a given time window [a i , b i ] within which the node may be visited. The SPPTW consists in finding the least cost route between two nodes ′ s ′ and ′ t ′ in the graph while visiting each node in the specified time window, it can be formulated as follows:

(SP P T W )                    min (i,j)∈A c ij x ij j∈V x sj = j∈V x jt = 1 j∈V x ij -j∈V x ji = 0, ∀i ∈ N \ {s, t} x ij ∈ {0, 1}, ∀(i, j) ∈ A x ij (T i + t ij -T j ) ≤ 0, ∀(i, j) ∈ A T j ∈ [a j , b j ], ∀j ∈ N T j ∈ N, ∀j ∈ N
where x ij , (i, j) ∈ A is the flow binary variable and T i , i ∈ N is the visiting time of node i.

Dynamic programming algorithm

We associate to each l-th path X l sj from the source s to the node j, a (time, cost) label denoted by E l j = (T l j , C l j ), corresponding respectively to the duration and the cost of X l sj . This label is computed iteratively along the path X l sj = (v 0 , v 1 , ..., v q ), where v 0 = s and v q = j as follows:

T v 0 = a v 0 , C v 0 = 0 T v i = max{a v i , T v i-1 +tv i-1 v i }, i = 1, ..., q C v i = C v i-1 + c v i-1 v i , i = 1, ..., q.
A label

E l j is feasible iff T v i ≤ b v i , ∀i = 1, .
.., q. Dominance rules are used to identify and discard useless labels. Definition 1. Let be X 1 sj and X 2 sj two different feasible paths from s to j and

E 1 j = (T 1 j , C 1 j ), E 2 j = (T 2 j , C 2 j ) the associated labels respectively. Then, E 1 j dominates E 2 j iff (T 2 j , C 2 j ) ≥ (T 1 j , C 1 j ) and (T 2 j , C 2 j ) = (T 1 j , C 1 j ).
If E 1 j is not dominated by another label, it is called efficient.. The complete set of efficient labels characterize the Pareto optimal frontier (Figure 1). We denote by E j = {E 1 j , E 2 j , ..., E p j } the set of efficient labels computed, and E p+1 j a new label at node j ∈ N . We present on Algorithm 1, the main steps of the dominance procedure applied on E p+1 j .

Algorithm 1 Dominance operations applied on a new label in a dynamic programming procedure next ← true, r ← 1.

repeat if E p+1 j = E r j then E j = E j ∪ {E p+1 j }, next = false. if E p+1 j > E r j then next ← false. if E p+1 j < E r j then E j = E j \ {E r j }, r ← r + 1. else r ← r + 1 until ((next = false) or (r > p)) if next then E j = E j ∪ {E p+1 j }
The treatment of a great number of labels increases the computation time of the dominance procedure. We present in the next section an improvement method which the main goal is the reduction of the number of labels saved and the number of labels compared in the dominance procedure.

Dynamic programming with blocs

The principle of this method is the computation of spaces which contain only efficient labels. We associate to each node i a set of blocs, each bloc is defined by a triplet

B k i = (T k i(low) , T k i(upp) , C k i(upp)
) where T k i(low) corresponds to the lower bound on time, T k i(upp) define the upper bound on time, and

C k i(upp)
represents the upper bound on cost of the bloc indexed by k. We suppose that the blocs are not lower bounded by the cost.

Initialization

Let be i the treated node, the set of blocs is initialized by the bloc B 0 i = (a i , b i , +∞) (figure 2). The bloc B 0 i contains only feasible and efficient labels. Let be 3). Initially, none efficient label is saved at node i, E 1 i is then efficient. The new label E 1 i allows to characterize the new bloc 3), the bloc B 0 i is then updated

E 1 i = (C 1 i , T 1 i ) a new label computed at node i ∈ N (Figure
B 1 i = (T 1 i , b i , C 1 i ) (Figure
B 0 i = (a i , T 1 i , +∞).

Expansion and dominance

Proposition 1. Let be E 1 i , ..., E l i the set of efficient labels computed at node i and B i = {B 1 i , ..., B l i } the associated blocs.. We consider

E i = (C i , T i ) a new a i b i B 0 i = (a i , b i , +∞) E 1 i = (C 1 i , T 1 i ) C i C 1 i T 1 i T i
Figure 2: Initialization 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000 0000000

T 1 i a i b i B 0 i = (a i , T 1 i , +∞) B 1 i = (T 1 i , b i , C 1 i ) C i E 1 i = (T 1 i , C 1 i ) T i C 1 i
Figure 3: Label expansion label computed at node i ∈ N . If the label E i belongs to a bloc B p i , p = 0, ..., l (i.e.. T p i(low) < T i ≤ T p i(upp) and T i < C p i(upp) ), then:

1. E i is a feasible label.
2. E i is an efficient label.

The labels dominated by

E i are E p+1 i , ..., E q i which satisfy C i ≤ C p+1 i , ..., C i ≤ C q i respectively. proof 1.
1. Each bloc in B i ∪ B 0 i belongs to the initialized bloc B 0 i defined in section 4.1. All labels belonging to these blocs are then feasible. i , all efficient labels of node i have durations greater than T i , so the label E i is efficient.

(a) If

3. If E i ∈ B p i , p = 0, ..., l, then T i ≥ T p-1 i > T p-2 i ... > a i , we conclude that none label in the set {E 1 i , ..., E p-1 i
} can be dominated by E i . The labels which can be dominated by E i belong to the set {E p+1 i , ..., E l i }, more precisely, the labels dominated by E i are E p+1 i , ..., E q i where C i ≤ C p+1 i , ..., C i ≤ C q i respectively.

An illustrative example

Figure 4 presents an illustration of the dynamic programming with blocs on an example where l = 6. Let be E 7 i = (T 7 i , C 7 i ) a new label computed at node i. Suppose that E 7 i belongs to the bloc B 3 i , it is feasible and efficient. Labels dominated by E 7 i belong to the set {E 4 i , E 5 i , E 6 i } and satisfy C 7 i ≤ C q i , q ∈ {4, 5, 6}, i.e E 4 i . This label and the associated bloc B 4 i are discarded, the bloc

B 7 i = (T 7 i , T 5 i , C 7 i ) is saved and the bloc B 3 i is updated (B 3 i = (T 3 i , T 7 i , C 3 i )).
00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 00 We denote by B j = {B 0 j , E 1 j , ..., E l j } the set of blocs computed at node j ∈ N , and E p+1 j a new label at node j. We present on Algorithm 2, the main steps of the dominance procedure applied on E p+1 j in our new procedure.

11
B 0 i = (a i , T 1 i , +∞) B 4 i = (T 4 i , T 5 i , C 4 i ) B 6 i = (T 6 i , b i , C 6 i ) B 5 i = (T 5 i , T 6 i , C 5 i ) B 3 i = (T 3 i , T 4 i , C 3 i ) B 2 i = (T 2 i , T 3 i , C 2 i ) B 1 i = (T 1 i , T 2 i , C 1 i ) a i T 1 i T 2 i T 3 i T 4 i T 5 i T 6 i b i T i C i E 3 i = (C 3 i , T 3 i ) E 2 i = (C 2 i , T 2 i ) E 1 i = (C 1 i , T 1 i ) E 7 i = (C 7 i , T 7 i ) E 4 i = (C 4 i , T 4 i ) E 5 i = (C 5 i , T 5 i ) E 6 i = (C 6 i , T 6 i ) B 0 i = (a i , T 1 i , +∞) B 3 i = (T 3 i , T 4 i , C 3 i ) B 2 i = (T 2 i , T 3 i , C 2 i ) B 1 i = (T 1 i , T 2 i , C 1 i ) a i T 1 i T 2 i T 3 i T 5 i T 6 i b i T i T 7 i B 5 i = (T 5 i , T 6 i , C 5 i ) B 6 i = (T 6 i , b i , C 6 i ) B 7 i = (T 7 i , T 5 i , C 7 i ) C i E 3 i = (C 3 i , T 3 i ) E 2 i = (C 2 i , T 2 i ) E 1 i = (C 1 i , T 1 i ) E 7 i = (C 7 i , T 7 i ) E 5 i = (C 5 i , T 5 i ) E 6 i = (C 6 i , T 6 i )

Experimentations

The results reported on Table 1 concern the resolution of the first pricing problem (SPPTW) in a column generation process for the resolution of the vehicle routing problem with time windows on acyclic graphs. The capacity constraints are relaxed. We consider 4 Solomon test instances (C101 25 (25 costumers), C101 50 (50 costumers), C101 (100 costumers) and C1 2 1 (200 costumers)) and 40 randomly generated instances. 10 instances are generated Algorithm 2 Dominance operations applied on a new label in the dynamic programming with blocs procedure if E p+1 j belong to a bloc B s j ∈ B j then

E j = E j ∪ {E p+1 j }. if E p+1 j = E s j then the bloc B s j is updated (section 4.3). next ← true. r ← s + 1. repeat if C p+1 j <= E r j then E j = E j \ {E r j }, B j = B j \ {B r j }, r ← r + 1. else next = false. until ((next = false) or (r > p)) B j = B j ∪ {B p+1 j } (B p+1 j is the bloc computed throw E p+1 j (section 4.3)).
for each size: G 100 (100 costumers), G 120 (120 costumers), G 140 (140 costumers) and G 160 (160 costumers), all results reported for each randomly generated class size are average values over 10 test instances. All instances are solved by the label correcting algorithm 1. Dynamic programming with blocs (DP Blocs) decrease the number of labels treated by 79% compared to the label correcting algorithm (DP LC), this is due to the reduction of the space search of new efficient labels using the blocs. DP Blocs permits also to reducing efficiently the number of labels compared in the dominance procedure (more than 99%) compared to DP LC. These results are predictable, as showed in proposition 2, dominance procedure is applied on a restricted number of labels. These results show the efficiency of the DP Blocs which permits of decreasing the resolution time of the DP LC method by 88%.

Solomon

Conclusion

Traditionally used dominance rules in dynamic programming consists, after the expansion of a label, to verify if the new label is feasible and not dominated (in the opposite case, it is ignored) by comparing it with all efficient labels already computed. This procedure can be time consuming when the number of labels is large.

We propose in this paper an improving technique which objective is the reduction of the number of labels saved and the computation time of the dominance procedure. In the proposed approach feasibility and dominance tests are included in one, which consists to check if the treated label belongs to a bloc. In the case where the label belongs to a bloc, this method permits a direct access to dominated labels.

Experimentations on SPPTW instances reveal that this method reduces significantly the number of labels saved and the computation time of the dominance procedure.

Figure 1 :

 1 Figure 1: Pareto optimal frontier at node j

Figure 4 :

 4 Figure 4: Illustrative example of expansion and dominance procedures

Table 1 :

 1 Dynamic programming with blocs

			C101 25			C101 50
	instances	resolT	nbLab nbLabDom	resolT	nbLab	nbLabDom
	DP LC	0,1"	343	680	3,3"	1 673	2 695
	DP Blocs	0.04"	79	4	0,2"	204	15
	Solomon		C101			C1 2 1
	instances	resolT	nbLab nbLabDom	resolT	nbLab	nbLabDom
	DP LC	72,0"	5 941	10 510 2 526,0"	39 922	1 946 876
	DP Blocs	3,2"	491	68	87,0"	1 036	11 425
	Generated		G 100			G 120
	instances	resolT	nbLab nbLabDom	resolT	nbLab	nbLabDom
	DP LC	186,0" 200 586 21 700 018	426,0" 335 421	116 011 504
	DP Blocs	18,9"	14 446	98 912	24,1"	92 325	123 413
	Generated		G 140			G 160
	instances	resolT	nbLab nbLabDom	resolT	nbLab	nbLabDom
	DP LC	462,0" 454 729 21 220 354	138,0" 105 478	11 836 750
	DP Blocs	27,0"	17 366	94 812	17,1"	57 758	115 377
	resolT: resolution time.
	nbLab: number of labels treated (labels on which the dominance procedure is applied).
	nbLabDom: number of labels compared in the dominance procedure.


Nora Touati

Touati Nora

Lucas Létocart

Anass Nagih

Diversification and reoptimization procedures in column generation for the resolution of the acyclic vehicle routing problem with time windows

Keywords: Diversification, Reoptimization, Column generation, Vehicle routing problem with time windows .52 i0 i3 i7 i8

Diversification and reoptimization procedures in column generation for the resolution of the acyclic vehicle routing problem with time windows

Introduction

A well known strategy to reduce the number of iterations in CG, is to return to the Master Problem (MP) many negative marginal cost columns at each iteration (minimization problem case). Generally, this intensification typically overloads the MP, while the final optimal base contains a very restricted subset of all the generated columns. Diversification methods which consist of computing complementary columns at each iteration of CG, permit to decrease the number of generated columns by characterizing quickly a good approximation of the Pricing Problem's (PP) convex hull. This approach is then more efficient on the first iterations [START_REF] Touati | Solutions diversification in a CG scheme[END_REF][START_REF] Touati | Performances improvement of the column generation algorithm: Application on the vehicle routing problem with time windows[END_REF].

Another way to improve the CG algorithm is to solve efficiently the PPs, known to be very greedy on computing time in many applications. We interest for this purpose to reoptimization techniques, which consist of using some informations about the PP's resolution at an iteration k, for solving efficiently the PP at iteration k + 1 [START_REF] Touati | Reoptimization techniques in a CG scheme[END_REF][START_REF] Touati | Performances improvement of the column generation algorithm: Application on the vehicle routing problem with time windows[END_REF].

This work focuses on the combination of diversification and reoptimization techniques in CG. We apply here a diversification method on the first iterations, and a reoptimization procedure on the last ones. The relevance of this approach is demonstrated by computational experiments on Solomon and on randomly generated instances of the Acyclic VRPTW (AVRPTW).

Vehicle routing problem with time windows

The VRPTW ( [START_REF] Desrochers | A new optimization algorithm for the vehicle routing problem with time windows[END_REF], [START_REF] Desrosiers | Routing with time windows[END_REF], [START_REF] Fisher | Vehicle routing with time windows : Two optimization algorithms[END_REF]) involves the design of a set of minimum cost routes, originating and terminating at a central depot, and visiting each costumer once, for a fleet of identical vehicles. The oriented network used by vehicles is defined by G = (V, A), where V is the set of nodes composed of the set of costumers N , a source 's' and a destination 't' representing exit and entrance to the depot respectively, and A the set of arcs connecting the customers and the depot. With each arc (i, j) ∈ A is associated a positive duration t ij and a cost c ij . Each vehicle has a given capacity Q and each customer i has a demand d i , i ∈ N . For each customer i, the start of the service must be within a given time window [a i , b i ]. A vehicle can wait in case of early arrival. The VRPTW can be formulated as a set partitioning problem, its linear relaxation is called the master problem and given as follows :

(M P V RP T W ) min r∈R c r λ r , r∈R δ ir λ r = 1, ∀i ∈ N , λ r ≥ 0, ∀r ∈ R
, where R is the set of all feasible paths, i.e. paths going from 's' to 't' and satisfying capacity and time windows restrictions, c r is the cost of path r and δ ir is 1 if path r visits node i and 0 otherwise. The PP is equivalent to a Shortest Path Problem with Time Windows and Capacity Constraints (SPPTWCC) and formulated as follows:

(SP P T W )                        min (i,j)∈A (c ij -π i )x ij j∈N x sj = j∈N x jt = 1 j∈N x ij -j∈N x ji = 0, ∀i ∈ N j∈N d j i∈N ∪{s} x ij ≤ Q x ij ∈ {0, 1}, ∀(i, j) ∈ A x ij (T i + t ij -T j ) ≤ 0, ∀(i, j) ∈ A T j ∈ [a j , b j ], ∀j ∈ N T j ∈ N, ∀j ∈ N
where x ij , (i, j) ∈ A are binary flow variables taking the value 1 when arc (i, j) is used by a vehicle, T i , i ∈ N are time variables associated with the starting of the service at node i, and π i is the dual variable associated with the i-th constraint of the MP. The standard approach to solve the SPPTW in practice is based on the dynamic programming method and has a pseudo-polynomial complexity. The principle is to associate with each possible partial path a label indicating the cumulated cost, time and demand, and to eliminate labels with the help of dominance rules. We consider here a label correcting algorithm, where nodes are repeatedly treated and their labels extended [START_REF] Desrochers | An algorithm for the schortest path problem with ressource constraintes[END_REF]. Solving the SPPTW with a dynamic programming approach is closely related to solving multi-objective shortest path problems, the aim is to generate efficient paths (i.e., Pareto optimal paths).

The main goal of our study is to evaluate the impact of two accelerating techniques of the column generation algorithm. Our study focuses on a particular version of the vehicle routing problem with time windows, based on acyclic networks, nevertheless many routing problems in practice are modeled by acyclic networks, particularly when the width of time windows is lower than the intertask time. It is the case in rail transportation problems for example.

Vehicle routing problem with time windows test instances

We consider in our experimentations two different types of AVRPTW instances. Some Solomon test problems from which we can extract acyclic networks (C101 25 (25 nodes), C101 50 (50 nodes), C101 (100 nodes), C1 2 1 (200 nodes)), and 10 AVRPTW instances uniformly randomly generated for each size: G 100 (100 nodes), G 120 (120 nodes), G 140 (140 nodes) and G 160 (160 nodes).

Diversification in a column generation algorithm

It is well known that the intensification allows to decrease the number of iteration of CG. Experimentations on AVRPTW instances described in section 2, reveal two interesting phenomena. The average percentage of suboptimal columns in the final optimal base is higher than 93% and 99% of the generated columns do not belong to the final optimal base.

To avoid the fast and useless increase of MPs size, we classically restrict to the generation of the k best solutions, k is a parameter expressed in % in what follows. We call this procedure k% Intensified Column Generation (k% ICG). Intensified Column Generation (ICG) correspond to k = 100. Experimentations on the AVRPTW [START_REF] Touati | Solutions diversification in a CG scheme[END_REF] showed that in average, the number of columns increases when k increases, whereas the number of iterations and the global resolution time decrease. Even if the add of columns with good reduced cost decreases the global number of generated columns, it increases considerably the number of iterations and is not enough to improve the global computing time. Instead of select columns based on the reduced cost criterion, we interest in what follows to columns structure based criterion.

Suboptimal columns characteristics

The previous experimental results show that suboptimal columns of best reduced cost are neighbor compared to the optimal column obtained, i.e. they covered almost the same tasks. We expose bellow the 20-th columns of best reduced cost, generated for a Solomon test instance of 25 costumers. The first column represent the path's cost and the others, the sequence of path's nodes.

In the PP's feasible solutions space, these columns are close compared to the PP's optimal solution, so, they can not improve the MP's feasible convex hull approximation. In the dual space, cuts corresponding to close columns are redundant for the current dual function approximation.

Thus, the addition of very close columns can be useless. The ICG needs to generate less columns, but good ones, for improving significantly the MP's convex hull approximation (and the dual function approximation). For this purpose, instead of selecting columns according to the reduced cost criterion, we interest in column's structure criterion. We look for complementary columns (i.e. columns don't having joint nodes) characteristics.

It is known that the add of complementary columns (resp. orthogonal cuts) to the MP permits to widen its convex hull (resp. the dual function) approximation [START_REF] Desrochers | A reoptimization algorithm for the shortest path problem with time windows[END_REF][START_REF] Vanderbeck | Decomposition and Column Generation for Integer Programs[END_REF]. It is showed in [START_REF] Touati | Solutions diversification in a CG scheme[END_REF] that this approach is more efficient on the first iterations when it is necessary to compute quickly a good primal (or dual) solution. On the last iterations of CG, the MP's convex hull is well approached and it is preferable in this case to generate columns to complete the current basis in order to obtain an optimal solution.

Diversification by resolution procedure and experimentations

Diversification procedures consist of inserting at each iteration of CG, a set of complementary columns to the MP. We consider here the CG with Diversification by Resolution (CGDR) [START_REF] Touati | Solutions diversification in a CG scheme[END_REF], which consists at each iteration of CG to solve a sequence of SPPTWCC, each one providing a path. After solving the l-th problem for a set of nodes N l , the nodes S l belonging to the optimal solution obtained are removed from N l to form a set of nodes N l+1 = N l \ S l which will be used for the (l + 1)-th problem.

We applied this procedure on AVRPTW instances described in section 2. Diversification by resolution is applied on the first iterations, when v(M P k )-v(M P k-1 ) 1 shows that diversification permits to decrease the number of computed columns without increasing considerably the number of iterations (compared to 20% ICG and 50% ICG). We generate with CGDR in average 61%, 69% and 79% of columns generated by ICG, 50% ICG and 20% ICG procedures respectively, this decreases MPs size and their resolution time. Diversification reduces the ICG resolution time by 39%.

v(M P k ) ≥ ε,

A reoptimization technique in CG

Consider the following problems:

(P ) : {max cx, x ∈ X ⊂ N n } and (P ′ ) : {max c ′ x, x ∈ X ′ ⊂ N n },
where c, c ′ ∈ R n , n ∈ N. Reoptimization consists of studying possibilities of an efficient resolution of (P ′ ) using some informations about the resolution of (P ).

In the CG process, only some PP's objective coefficients change, due to the update of MP's dual solutions. So, reoptimization in this case can be used. The principle is to retain some informations from the PP's resolution at an iteration k, to be used at iteration k + 1, for accelerating the PP's resolution.

We denote by E k j the set of efficient labels of node j ∈ V at iteration k of CG. The set E k j permits to build a Pareto frontier P k j representing a lower bound on the cost of all feasible labels at node j. Master problem's resolution at iteration k + 1, gives new dual variables values. Labels in E k j are still feasible at this iteration but their cost change, they define a primal function representing an upper bound on costs of labels in E k+1 j as function of time. Then, the Pareto optimal frontier computed at iteration k allows to define a primal function at iteration k + 1. To compute the set E k+1 j , it is enough to calculate all labels dominating those in E k j .

Labels computed at iteration k + 1 belong to the space between the primal function and P k+1 j , called treatment zones. We reduce then the space search of efficient labels, the number of labels computed and the dominance procedure time. The smaller are the treatment zones, the smaller is the space search of new efficient labels and more reoptimization is efficient.

It is frequently observed in a CG process, that MP's values evolve slightly on the last iterations, this phenomenon is called the long tail effect, so, the dual variables are also close, and the change of arcs costs is weak. The treatment zones on these iterations can be small and reoptimization is more efficient.

We tested this approach on AVRPTW instances. ICG is used until the inequality 2 show that in average, GCReopt decrease the number of computed labels by 53% compared to ICG. This is due to the use at each iteration of efficient labels computed at the previous iteration, allowing the construction of primal functions which restrict the space search of new efficient labels. Pricing problems resolution time is then significantly decreased. The average gain of computing time obtained on the global resolution time of the ICG algorithm is 36%. We observe that the gain obtained by the GCReopt on the ICG global computing time is not as important as the gain on the number of treated labels, this is due to the overcost of primal functions construction.

v(M P ) k -v(M P ) k-1 v(M P ) k ≤ ε is satisfied, next,

Diversification and reoptimization in a CG algorithm

In the CG algorithm we can distinguish two important phases. The first one, is the begin of the process, where the MP is poor in informations. Many methods in literature propose techniques to enrich the MP (start the CG with a pool of columns generated heuristically or with a good dual solution ...). We interest here on using complementary solutions to provide good informations for the MP by characterizing efficiently a good approximation of its convex hull. The second one is the end of the process, where we frequently observe a weak evolution of MPs values. On the last iterations the dual variables values are neighbor, so, reoptimization methods can be efficient to reduce the resolution time of the PP.

We showed that the ICG is improved by using diversification or reoptimization techniques. We aim to study the impact of their combination on CG performances. We use diversification on the first iterations until the inequality v(M P ) k -v(M P ) k-1 v(M P ) k >= ε is satisfied and reoptimization on the other iterations. We call this procedure Column Generation with Diversification and Reoptimization (CGDReopt). Table 3 presents a comparative study between ICG, CGDR, CGReopt and CGDReopt procedures, by experimentations on the test instances described in section 2. The ICG procedure computing time is improved on one hand by using diversification on the first iterations where the gain obtained is in average 48% and on another hand by using reoptimization on the last iterations, this decrease the ICG resolution time by 34% in average. The gain on the global computing time is related to the number of iterations made with diversification and reoptimization respectively.

As expected, the use of both these two methods decreases significantly the resolution time of the ICG procedure, the average gain obtained on computing time is respectively.

For our experimentations, the compromise between the cost of treatment of the accelerating procedure (the computation of complementary columns for CGDR and the construction of primal functions in CGReopt) and the gain on the global resolution time (gain on the number of generated column in CGDR and the gain on the number of treated labels for CGReopt) is better for the CGDR procedure than CGReopt.

Conclusion

Diversification methods give a very interesting gain on the global number of generated columns in GC by building a good approximation of PP feasible domain. It reduces the MP's size and the global resolution time. These approaches are more interesting when used on the first iterations where the PP's convex hull approximation is poor.

Reoptimization approaches are used successfully to improve PP's resolution time in CG. We present a reoptimization algorithm in CG for the resolution of AVRPTW. The principle is to use efficient labels computed at an iteration k to construct a primal function for reducing the space search of new efficient labels at iteration k + 1. This reduces the number of treated labels and dominance time.

We combine in this work these two alternatives in CG by applying a diversification method on the first iterations and a reoptimization procedure on the last ones. Experimental results on AVRPTW instances show that the gain on computing time obtained by this combination is much more important than that obtained by the two approaches used independently.

Table 1 :

 1 where v(M P k ) is the MP's value at iteration k. The intensification is used on the other iterations. All results reported in table 3 for each randomly generated class size are average values over 10 test instances. Generation of complementary solutions Table

	Solomon's		C101 25			C101 50			C101				C1 2 1
	instances	nbIt	nbC	GRT	nD	nbI	nbC	GRT	nD	nbIt	nbC	GRT	nD	nbIt	nbC	GRT	nD
	20% ICG	51	153	4,3"	0	127	647	73,1"	0	170	2 008	15,2'	0	286	10 768	438,6'	0
	50% ICG	30	201	2,8"	0	61	823	32,7"	0	105	2 154	9,6'	0	250	11 020	197,1'	0
	ICG	24	246	2,1"	0	55	1 003	31,7"	0	63	2 724	5,0'	0	163	11 460	150,9'	0
	CGDR	35	251	3,0"	3	83	956	51,8"	4	75	2 646	2,1'	2	114	10 141	102,2'	5
	Generated instances	nbIt	G 100 nbC	GRT	nD	nbIt	G 120 nbC	GRT	nD	nbIt	G 140 nbC	GRT	nD	nbIt	G 160 nbC	GRT	nD
	20% ICG	69	6 263	10,3'	0	99	7 864	32,3'	0	97	10 169	43,5'	0	72	10 409	18,0'	0
	50% ICG	36	11 406	9,8'	0	48	13 663	31,2'	0	49	18 800	40,1'	0	42	17 980	20,0'	0
	ICG	27	16 960	12,5'	0	26	19 013	39,1'	0	25	25 259	35,5'	0	27	23 864	20,4'	0
	CGDR	40	2 589	5,2'	22	50	2 579	17,5'	28	55	5 370	15,4'	37	48	1 298	6,8'	31
	nbIt: global number of iterations.
	nbC: number of generated columns.
	GRT: global resolution time ((') minutes and (") seconds).
	nD: number of iterations with diversification.
	ε = 0.001.

  the reoptimization technique is used. We call this procedure CG with Reoptimization (GCReopt).

	Solomon's		C101 25			C101 50			C101			C1 2 1
	instances	GRT	nbLb	ItR/ItG	GRT	nbLb	ItR/ItG	GRT	nbLb	ItR/ItG	GRT	nbLb	ItR/ItG
	ICG	2,1"	9 112	0/24	31,7"	86 073	0/55	5,0'	384 810	0/63	150,9'	4 262 215	0/163
	GCReoptP	1,2"	3 968	18/24	29,0"	18 178	48/55	2,1'	118 292	52/63	76,8'	1 148 708	114/163
	Generated instances	GRT	G 100 nbLb	ItR/ItG	GRT	G 120 nbLb	ItR/ItG	GRT	G 140 nbLb	ItR/ItG	GRT	G 160 nbLb	ItR/ItG
	ICG	12,5'	446 080	0/27	39,1'	5 841 413	0/26	35,5'	454 729	0/25	20,4'	475 701	0/27
	GCReoptP	8,9'	315 241	8/27	21,2'	611 876	9/26	21,1'	405 266	7/25	17,3'	423 800	4/27
	GRT: global resolution time ((') minutes and (") seconds).
	nbLb: number of treated labels.
	ItG: global number of iterations.
	ItR: number of iterations with reoptimization.
	ε = 0.001.

Table 2 :

 2 Reoptimization in a CG algorithmExperimental results reported on table

Table 3 :

 3 Combination of diversification and reoptimization in CG