Bernard Bonaïti

Valérie Bonadona

Hervé Perdry

Nadine Andrieu

Catherine Bonaïti-Pellié
email: catherine.bonaiti@inserm.fr

Estimating penetrance from multiple case families with predisposing mutations: Extension of the "genotype-restricted likelihood" (GRL) method

Keywords: penetrance, bias, pleiotropy, parent-of-origin, families

Some diseases are due to germline mutations in predisposing genes, such as cancer family syndromes.

Precise estimation of the age-specific cumulative risk (penetrance) for mutation carriers is essential for defining prevention strategies. The genotype-restricted likelihood (GRL) method is aimed at estimating penetrance from multiple case families with such a mutation. In this paper, we proposed an extension of the GRL to account for multiple trait disease and to allow for a parent-of-origin effect.

Using simulations of pedigrees, we studied the properties of this method and the effect of departures from underlying hypotheses, misspecification of disease incidence in the general population or misspecification of the index case, and penetrance heterogeneity. In contrast with the previous version of the GRL, accounting for multiple trait disease allowed unbiased estimation of penetrance. We also showed that accounting for a parent-of-origin effect allowed a powerful test for detecting this effect.

We found that the GRL method was robust to misspecification of disease incidence in the population, but that misspecification of the index case induced a bias in some situations for which we proposed efficient corrections. When ignoring heterogeneity, the penetrance estimate was biased toward that of the highest risk individuals. A homogeneity test performed by stratifying the families according to the number of affected members was shown to have low power and seems useless for detecting such heterogeneity. These extensions are essential to better estimate the risk of diseases and to provide valid recommendations for the management of patients.

INTRODUCTION

Some diseases with variable age of onset are due to the presence of predisposing gene mutations, such as mismatch repair (MMR) genes in hereditary nonpolyposis colorectal cancer (HNPCC) or BRCA1 and BRCA2 in breast-ovarian cancer syndrome. These genes may be responsible for hereditary forms of these diseases. Precise estimation of the age-specific cumulative risk (penetrance function) for mutation carriers is essential for defining prevention strategies.

Families in which such mutations have been identified can contribute to estimate these risks, as long as adjustment is made for these families generally ascertained because of several affected members [START_REF] Carayol | Hereditary non-polyposis colorectal cancer: current risks of colorectal cancer largely overestimated[END_REF] . Such families are usually referred by physicians to genetic counsellors who propose genetic testing when specific criteria are fulfilled. For hereditary cancers, for example, most criteria used for recommending genetic testing are based on familial aggregation of specific cancers [START_REF] Eisinger | Identification et prise en charge des prédispositions héréditaires aux cancers du sein et de l'ovaire[END_REF][START_REF] Lindor | Recommendations for the care of individuals with an inherited predisposition to Lynch syndrome: a systematic review[END_REF] . When a mutation is identified in an affected member (defined as index case), genetic testing is proposed to close relatives who, if carriers, will be offered intensive surveillance, which will improve the prognosis, or prophylactic surgery when possible.

An ascertainment-adjusted method, based on maximum likelihood, has been proposed for estimating the age-specific cumulative risk (penetrance) of a given disease associated with a deleterious mutation from families in which such a mutation has been identified [START_REF] Carayol | Estimating penetrance from family data using a retrospective likelihood when ascertainment depends on genotype and age of onset[END_REF] . This method, called the "genotype-restricted likelihood" (GRL) method, provides unbiased penetrance estimates whatever criteria used for ascertainment of families and without having to model the ascertainment process. The GRL method corrects for the bias due to the selection on carrier genotype of the index case since only families with an identified carrier individual are informative for penetrance estimation. This method is especially appropriate for hereditary predispositions to common diseases when numerous and complex familial criteria involving several affected relatives are used to recommend genetic testing. It has been shown to be independent of selection criteria, in particular on the number of affected individuals and on the age at diagnosis of affected family members [START_REF] Carayol | Estimating penetrance from family data using a retrospective likelihood when ascertainment depends on genotype and age of onset[END_REF] .

Beside the numerous advantages mentioned above, the GRL method relies on assumptions which may not be fulfilled in some situations. In order to evaluate its robustness to a departure from these hypotheses, we studied the sensitivity of the GRL in various situations such as disease frequency of the general population over-or under-estimated or a genetic heterogeneity not taken into account.

Additionally, the previous GRL version allows estimating the penetrance of only one trait once at a time and this might bias penetrance estimates in family syndromes where several different traits may occur (pleiotropy), like different tumor localizations. Therefore, we proposed in this paper to extend the method to account for pleiotropy. We also extended the method to allow for a parent-of-origin effect, i.e. penetrance functions differing according to the gender of the parent who transmitted the deleterious mutation. This effect has been described in some diseases [START_REF] Shete | Genetic imprinting analysis for alcoholism genes using variance components approach[END_REF][START_REF] Green | Impact of gender and parent of origin on the phenotypic expression of hereditary nonpolyposis colorectal cancer in a large Newfoundland kindred with a common MSH2 mutation[END_REF][START_REF] Gorlova | Imprinting detection by extending a regression-based QTL analysis method[END_REF][START_REF] Klutz | A parent-of-origin effect in two families with retinoblastoma is associated with a distinct splice mutation in the RB1 gene[END_REF] but has not been yet accounted for in penetrance estimation. These extensions are essential to better estimate the risk of diseases and to provide valid recommendations for the management of patients.

METHODS

The GRL method

The GRL [START_REF] Carayol | Estimating penetrance from family data using a retrospective likelihood when ascertainment depends on genotype and age of onset[END_REF][START_REF] Alarcon | Estimating cancer risk in HNPCC by the GRL method[END_REF] uses a retrospective likelihood conditioned on the phenotypes of all family members and on the genotype of the index case. It estimates penetrance parameters for mutation carriers by maximizing the probability of observed genotypes (G) of family members who have been tested for the mutation found in the index case, conditional on observed phenotypes (P) and on index case being a carrier (I). Due to the fact that the index case is always tested, the conditional probability may be written as: Pr(G/P,I) = Pr(G,P) / Pr(P,I) Following the cure model of De Angelis 10 , we considered that a proportion κ of individuals will never be affected and a Weibull model for penetrance for the others. The cumulative risk by age t is:

( ) [ ] { } l t exp ) ( ) , , ; t ( F ) t T Pr( l l l l l α λ κ λ α κ - - - = = < 1 1
where l is the genotype for the mutation (l = 1, 2, 3 for AA, Aa and aa genotypes respectively, A being the mutated allele), α l and λ l are the Weibull shape and scale parameters, κ l the probability of never being affected given l, and T is the age at disease occurrence.

Let y ij =(t ij , δ ij , g ij ) the set of observations on the j th member of the i th family, t ij the age at onset of the disease or the age at censoring (earliest date among dates of prophylactic surgery, death or last news), δ ij the indicator of occurrence of the disease before the age at censoring and g ij the observed genotype which is coded as 0 when unknown. Let F(t; θ l ) the cumulative risk by age t for the l th genotype and h(t; θ l ) the corresponding hazard functions with θ l =(α l , λ l , κ l ). The probability of the set of observations (y ij ) given the latent genotype (u ij ) is:

( ) ( ) ( ) ( ) Pr 1 , , , ij ij ij ij l ij l ij y u l F t h t g l δ θ θ ϕ ⎡ ⎤ ⎡ ⎤ = = - × × ⎣ ⎦ ⎣ ⎦
where φ(g ij ,l) is the probability of observed genotype knowing the latent genotype : φ(g ij ,l) = 1 if g ij =0 or l and 0 otherwise.

Finally it is possible to calculate the probability Pr(G,P) using

( )

Pr ij ij y u l = with the Elston- Stewart algorithm [START_REF] Elston | A general model for the genetic analysis of pedigree data[END_REF][START_REF] Thompson | Pedigree Analysis in Human Pedigrees[END_REF] . The denominator, Pr(P,I), is computed in the same way with g ij being known only for the index individual.

The maximization of the conditional likelihood Pr(G/P,I) was performed by the L-BFGS-B algorithm [START_REF] Zhu | L-BFGS-B: Algortithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization[END_REF] . When analysing a real family sample, confidence intervals of penetrance estimates are obtained by bootstraps.

The method assumes that the mutated allele frequency in the general population is known as well as the penetrance function in non carriers which is set equal to the incidence in the general population.

Extensions of the GRL

Multiple trait phenotype

The contribution to the likelihood of each individual was modified to simultaneously take into account the phenotype of a variable number of possible traits, under the hypothesis that, given genotype, the occurrence of a disease does not modify the risk of developing subsequently other diseases. If t ijk is the age at onset or the age at censoring of the k th disease and δ ijk the indicator of the occurrence of the k th disease before censoring time on the j th member of the i th family, y ij =(t ij1 t ijk .. , δ ij1

.. δ ijk .., g ij ) with probability :

[ ] [ ] { } ) l , g ( ) , t ( h ) , t ( F ) l u / y Pr( ij k kl ijk kl ijk ij ij ijk ϕ θ θ × × - = = ∏ ∂ 1

Parent-of-origin effect

To take into account a parent-of-origin effect, we modified the likelihood by splitting the heterozygote genotype in two different genotypes according to the paternal or maternal origin of the mutated allele. Four genotypes were considered: AA, Aa, aA and aa, where Aa and aA are the ordered heterozygous genotypes in which the first allele is transmitted by the father and the second one by the mother. The matrix of genotype probabilities from one generation to another was modified accordingly and three penetrance functions instead of two were considered for gene carriers.

Simulations of pedigrees

To study the statistical properties of the method, and in particular its robustness to departures from underlying hypotheses, as well as the interest of the extensions implemented in the method, we simulated four generation families with a fixed number of relatives: starting from the index case, we simulated a couple of his(her) parents and two couples of grand-parents. The grand-parents and parents of the index case (generation 1 and 2) had respectively two and three children. Each of them had two children as well as all their offspring until the fourth generation (figure 1). Genotypes of family members were randomly generated according to the mutated allele frequency in ancestors and spouses and to Mendel's laws for offspring except the index case who had at least one mutated allele.

To simulate realistic situations, we chose to compute parameters to fit the data of a large national French survey of 537 families with Lynch syndrome [START_REF] Bonadona | Cancer risks associated with germline mutations in MLH1, MSH2 and MSH6 genes in Lynch syndrome: results from the large nationwide French study ERISCAM[END_REF] . Age at last news (or age at death) for each individual was assumed to be normally distributed with 63, 60, 50, 35 years for means and 17, 16, 13 and 11 years for standard deviations respectively for generations 1 to 4. For each family member (including the index case), age at disease occurrence was randomly generated given penetrance function according to genotype. The individual was set as affected or not by comparing the age at disease occurrence and the age at last news. Finally, a family was selected if the index case was affected and if the family fulfilled selection criteria on the minimal number of affected (MNA) individuals.

We considered a dominant genetic model for disease transmission with equal penetrance functions for AA and Aa genotypes which has been often observed in cancer predisposition so far. The A allele was set to 0.001, so that most of the mutation carriers have the Aa genotype. The sets of penetrance functions used for the simulations are shown in table 1. Three sets of penetrance fonctions (high, medium and low risks) were used in the case when only one trait may be observed. To simulate realistic situations for the several different traits exercise, parameters for the simulation of families were chosen to fit the risks published for colorectal, endometrial, ovarian, and urologic tract cancers in Lynch syndrome [START_REF] Alarcon | Estimating cancer risk in HNPCC by the GRL method[END_REF][START_REF] Dunlop | Cancer risk associated with germline DNA mismatch repair gene mutations[END_REF][START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF][START_REF] Jenkins | Cancer risks for mismatch repair gene mutation carriers: a population-based early onset case-family study[END_REF][START_REF] Watson | The risk of extra-colonic, extra-endometrial cancer in the Lynch syndrome[END_REF] .

As only some family members are usually tested in families, we had also to simulate when genotypes are known or not, with probabilities fitting to realistic situations in cancer genetics. That is, we fixed some individuals as known or unknown for their genotype and generated the availability on the other family members' genotypes using a probability π, subsequently referred as the genotyping rate, of being tested as follows:

-The index case and his(her) parents are always tested -Grand-parents and spouses are never tested, as well as the sib of the non carrier parent and his (her) offspring -The sib of the carrier parent is tested with probability π and, if non carrier, his(her) offspring is not tested -The offspring of family members found to be carriers is tested with probability π.

Properties of the GRL

Estimates and standard errors of penetrance at different ages between 20 and 80 years were obtained from their average on 200 replicates in different situations according to penetrance among carriers and non carriers, genotyping rate π (0.1, 0.3, 0.5, 0.8, 1.0), number of families in the sample (100, 200, 500), and MNA (2, 3 or 4) individuals for a family to be ascertained. Maximum likelihood estimation should be asymptotically unbiased but some bias may exist with small sample size. It was measured by the average difference between estimations and the true value.

Robustness to departures from underlying hypotheses was investigated by invalidating these hypotheses and by evaluating the bias induced on penetrance estimates when analysing the simulated data. We considered two sources of error when using the GRL method: 1) misspecification of the disease incidence in the general population; 2) misspecification of the index case (cf. infra). In addition, the GRL assumes that the probability of phenotypes is the same for all mutation carriers, i.e. no heterogeneity in penetrance. However, heterogeneity of penetrance may exist because of modifier factors and one may expect that families with higher penetrance are likely to have more affected individuals than families with lower penetrance. We considered the possibility of penetrance heterogeneity and studied estimates obtained when ignoring this heterogeneity. We also studied the possibility of detecting heterogeneity using the number of affected individuals in the family as a surrogate. Homogeneity tests were performed using likelihood ratio tests that compare the maximum likelihood value L 1 obtained assuming the same penetrance in the different subgroups to the maximum likelihood value L 2 obtained when allowing for a difference among subgroups. The statistic used, -2 Ln [L 1 /L 2 ], follows a Chi square distribution with a number of degrees of freedom equal to the difference between the number of parameters estimated in L 1 and in L 2 respectively.

The interest of the extension to a possible parent-of-origin effect was evaluated by simulating such an effect in two scenarios i) strong effect with large difference in penetrance, ii) small effect with moderate difference in penetrance according to the parental origin of the mutation. In these two situations, we estimated the power to detect the effect using a likelihood ratio test with a 0.05 type I error. The interest of the extension to a multiple trait phenotype was evaluated by comparing the results of the joint analysis to those obtained by performing separate analyses for each trait, where individuals are considered as unaffected when affected by another trait than the one under study.

RESULTS

Bias and efficiency of the GRL

The bias on penetrance estimates due to limited sample size is shown in table 2 for various situations. The bias is most often positive (i.e. penetrance is overestimated) but very small. Only in extreme situations when the sample size is small (i.e. 100 families), few individuals are tested in the family (e.g. π=0.1) and penetrance is low, the bias may be non negligible compared to the penetrance. This bias increases with age due to the decrease in genotyped individuals among older family 9 members. There is no variation in bias according to the MNA individuals required for the selection of pedigrees (data not shown).

The efficiency, measured by the standard error of the estimates, also depends on the age, the sample size, the proportion of known genotypes and the true risk values. The results are given in table 2 and show small standard errors for most of estimates. For example, a sample of 200 families with a low genotyping rate (π=0.3) would provide a standard error of 0.07 by age 70 years for a disease with medium risk. The standard errors do not vary with the MNA individuals required for the selection of pedigrees (data not shown).

Robustness to departures from underlying hypotheses

Misspecification of the disease incidence in the general population

The incidence of common diseases in the general population is known when regional or national registries provide accurate rates. However, these registries are usually available for recent time periods only: e.g. the first cancer registry was established in France in 1978. As cancer incidence has varied with time, the incidence specified in the analysis may not be valid for the past generations. In many other diseases for which there are no registries, the incidence may be estimated with some error. We estimated the bias induced by a large error, i.e. by dividing or multiplying the risks by 2 in the analysis, for various values of population incidence. The results are given in table 3. The bias induced by an error on population incidence is rather small and less than 0.05 in a majority of situations. It barely varies with π and the MNA in the families (data not shown).

Misspecification of the index case

The index case is defined as the first individual who was tested in the family. Index cases are often incident cases who asked for genetic counselling because one or several relative(s) are also affected.

However, when there are only prevalent cases in a family, the geneticist has to choose, among the affected members, the person who will first be tested and will be considered as the index case. In practice, the geneticist chooses the one with the highest prior probability of being a carrier, in general the youngest affected one. Such a choice is expected to induce a bias on penetrance estimated as explained on a specific example given in the appendix. To investigate this potential bias, we defined as the index case the youngest affected one in the analyses. The results are shown in table 4 only for MNA=3 or 4. Indeed, no effect was observed for MNA=2 since, in the simulated data, the index case was often the youngest affected one. In general the bias is small, but may be substantial for high and medium penetrance values, with an underestimation of the penetrance for young ages and an overestimation for old ones.

To possibly correct for this bias, we proposed in table 4 two strategies of analysis that free to know the "true" index: 1) indicate as the index case in the analysis an individual chosen at random among affected family members, 2) condition the likelihood on the event that at least one affected individual is a mutation carrier as proposed by Quehenberger et al. [START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF] , with two options for affected individuals participating to the conditioning i) all affected individuals, tested or not; ii) only tested affected individuals. The three methods most often provide a decrease in bias compared with choosing the youngest one; globally, the smallest bias is obtained when choosing the index at random.

Regarding the second strategy, a smaller bias was obtained when using the second option, i.e. only tested affected individuals participate to the conditioning.

Penetrance heterogeneity

We considered the case where the population of carriers would be a mixture (50% each) of low risk (penetrance by age 40, 70 and end of life: .05, .20 and .22 respectively) and high risk families (penetrance by age 40, 70 and end of life: .20, .80 and .90 respectively). We analyzed the data under the assumption of homogeneity, and also by stratifying on the number of affected individuals (less than 3 vs. at least three, and less than 4 vs. at least 4) in the family and performed a homogeneity test at the 5% level. The power was estimated as the proportion of 400 replicates with a significant homogeneity test (table 5). As expected, the penetrance estimated when ignoring heterogeneity was toward the penetrance value of high risk families, due to selection on multiple case families. When we stratified the families according to the number of affected members, we obtained different risk estimates in the two strata, but the power for detecting heterogeneity was very low whatever the number of affected members used for stratification.

Extensions of the GRL

Multiple trait phenotype

We considered four different traits, with different values of penetrance between men and women for some traits in order to fit the risks of cancers described in Lynch syndrome as indicated in table 1.

The results are given in table 6. The penetrance estimates using single trait analyses are generally underestimated when compared with values in simulated data and the bias is stronger when the MNA individuals in the pedigrees is high, except for colorectal cancer in men. For instance, for endometrial cancer where the true penetrance by age 70 years was simulated at 0.28 and the penetrance estimated separately, i.e. by considering women affected by other cancers as unaffected, was respectively 0.27, 0.22 and 0.18 when families were ascertained on at least two, three, and four affected individuals, respectively. In contrast, the multiple trait analysis provided unbiased estimates.

Parent-of-origin effect

Table 7 shows the power to detect a parent-of-origin effect in samples of various sizes and variable genotyping rate π, for two sets of differences in penetrance according to the gender of the parent having transmitted the mutated allele: a large difference between penetrances by age 70 of respectively 0.30 and 0.60, and a small one of respectively 0.40 and 0.50. As for penetrance heterogeneity, the power was estimated as the proportion of replicates with a significant homogeneity test. The power to detect the parent-of-origin effect was found to be very high, even for samples of moderate size, when the difference was large, but decreased dramatically when the difference between penetrances decreased.

DISCUSSION

Since 1994, several genes have been identified the mutations of which are responsible for hereditary forms of cancer, in particular breast/ovarian cancer and colorectal cancer; guidelines for genetic testing and clinical management have been published and the detection of mutations is now routinely organized [START_REF] Schwartz | Proceedings of the international consensus conference on breast cancer risk, genetics, & risk management[END_REF][START_REF] Lynch | Hereditary colorectal cancer syndromes: molecular genetics, genetic counseling, diagnosis and management[END_REF] . In France, more than 2500 mutations in BRCA1 and BRCA2 were found in index cases and more than 7000 relatives had genetic testing between 2003 and 2007. For MMR genes, a mutation was found in more than 1000 index cases and nearly 3000 relatives were tested (http://www.e-cancer.fr/v1/fichiers/public/synthese_evolution_activite_2003_2007.pdf). Such familial data provide unique information on carrier risk, as long as adjustment for selection criteria of these families is adequately performed. These criteria are complex and have evolved with time.

Moreover, these criteria are proposed as guidelines for genetic counsellors and therefore are not purely applied. Therefore, a formal correction for these criteria is completely unfeasible and the GRL method appears particularly appropriate for estimating disease risks in carriers.

In this paper, we investigated the statistical properties of the GRL method and proposed some extensions.

As any maximum likelihood estimator, the GRL estimator is unbiased only asymptotically.

Indeed, we observed a small bias in penetrance estimate for samples of moderate size, particularly at old ages. This bias should be kept in mind when analyzing small samples of families. Regarding precision of the estimates, we showed that the standard errors in current situations were rather small, although retrospective likelihood methods are known to be poorly efficient [START_REF] Kraft | Bias and efficiency in family-based gene-characterization studies: conditional, prospective, retrospective, and joint likelihoods[END_REF][START_REF] Choi | Estimating disease risk associated with mutated genes in family-based designs[END_REF] .

We studied the robustness of the GRL to departures from underlying hypotheses. We found that the method was very robust to a misspecification of disease incidence in the general population, even for important errors. Such an error could be due to a cohort effect as observed in some cancers [START_REF] Belot | Cancer incidence and mortality in France over the period 1980-2005[END_REF] . Our results indicated that the method is expected to perform well when using an average disease incidence at an intermediate period between the periods of oldest and youngest generation of the families.

Regarding the specification of the index case, there is no ambiguity if the index is an incident case who asked for genetic counselling and the GRL should be used as it. If there is any doubt on the index identification in some families, when all affected family members are prevalent cases, our recommendation is to use, for these families only, either the random procedure or the conditioning on the genotype of all individuals affected by a disease under study [START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF] , provided that they had genetic testing proving their carrier status. In general, the random procedure provides the least biased estimates but Quehenberger's method [START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF] is a good alternative.

Regarding penetrance heterogeneity, we found that the estimated penetrance was close to that of highest risk individuals because of selection on multiple case families. A homogeneity test performed by stratifying the families according to the number of affected members was shown to have low power and appears useless for detecting such heterogeneity. However, one must keep in mind that nearly all families referred to genetic counselling so far are multiple case families and that the average penetrance estimated from a sample of such families may be the most appropriate one to predict disease occurrence and to recommend or not genetic testing in such context.

Another source of bias could be that we did not take into account a possible family-specific random effect due to low penetrance genes with multiplicative effects, such as those found by Antoniou et al. [START_REF] Antoniou | A comprehensive model for familial breast cancer incorporating BRCA1, BRCA2 and other genes[END_REF] in breast cancer, nor modifier genes such as those found to influence cancer risk in BRCA2 carriers [START_REF] Antoniou | Common breast cancerpredisposition alleles are associated with breast cancer risk in BRCA1 and BRCA2 mutation carriers[END_REF][START_REF] Antoniou | RAD51 135G-->C modifies breast cancer risk among BRCA2 mutation carriers: results from a combined analysis of 19 studies[END_REF] . Such an effect could affect major gene penetrance estimate [START_REF] Choi | Estimating disease risk associated with mutated genes in family-based designs[END_REF][START_REF] Carroll | Score tests for familial correlation in genotyped-proband designs[END_REF][START_REF] Li | Semiparametric estimation of major gene and family-specific random effects for age of onset[END_REF] . The GRL will have to be extended to take into account such effects when better identified.

The GRL method assumes that the frequency of deleterious mutations is low in the general population. A departure from this hypothesis would invalidate our approximation for risks in non carriers that are set to be equal to the risks in the general population. The sample size needed to estimate the risk in non carriers by the GRL is prohibitive and this approximation is necessary.

Therefore, we do not recommend using the GRL for mutations with frequency that would exceed 0.01 for a dominant predisposition, or 0.10 for a recessive one, in the general population.

Finally, we could evaluate the efficiency of the extensions of the GRL that we proposed, i.e.

possible multiple trait phenotype and parent-of-origin effect. Regarding the former, we found that analyzing separately each trait induced a bias which was corrected when using the multiple trait option. The bias is due to the fact that a single trait analysis uses an inaccurate ascertainment correction when conditioning on only one trait whereas several traits led to the ascertainment of families. Indeed, our results showed that the risks of CRC in men were unbiased when performing a single trait analysis, which is most probably due to the fact that CRC is almost the only tumour observed in carrier men. Therefore, this extension, in addition to sparing time by performing only one analysis, avoids bias, particularly for low risk traits. This is most probably the reason why we found a lower risk of endometrial cancer than other studies when we analyzed a sample of 36 French families

Figure 1 .

 1 Figure 1. Family structure used for the simulations (this is only an example as the sex of each individual is randomly generated)

Table 1 .

 1 Penetrance functions used for simulations

	phenotype Disease	gender Penetrance in mutation	Penetrance in non carriers
				carriers by age		by age
				40	50	70	40	50	70
	Single trait Low risk	all	0.01	0.02	0.05	0.0001 0.0006 0.005
		Medium risk all	0.10	0.19	0.40	0.001 0.006 0.05
		High risk	All	0.20	0.40	0.80	0.02	0.04	0.10
	Multiple	CRC	M	0.06	0.14	0.40	0.0009	0.013	0.032
	trait		F	0.04	0.10	0.30	0.0008	0.009	0.020
		END	F	0.05	0.11	0.28	0.0002 0.004 0.011
		OVA	F	0.02	0.04	0.07	0.001 0.006 0.010
		URE	M	0.005	0.009	0.02	0.0001	0.0005	0.0009
			F	0.005	0.009	0.02	0.0001	0.0003	0.0004
	M: male, F: female; CRC: colorectal cancer, END: endometrial cancer, OVA: ovarian cancer, URE:
	cancer of the urologic tract

Table 2 .

 2 Maximum likelihood bias (and standard errors) on penetrance estimate by ages 50 and 70 years, according to sample size and genotyping rate (π), for low, medium and high penetrance values (selection of families on at least two affected individuals, 200 replicates)

				Bias on penetrance estimate
					(standard error)
				By age 50		By age 70
	Penetrance π	Number of families	Number of families
			100	200	500	100	200	500
	Low	0.1	0.05	0.02	0.01	0.06	0.03	0.01
			(0.08)	(0.06)	(0.03)	(0.08)	(0.06)	(0.03)
		0.3	0.02	0.02	0.01	0.03	0.02	0.01
			(0.05)	(0.03)	(0.02)	(0.05)	(0.03)	(0.02)
		0.8	0.02	0.01	0.00	0.02	0.01	0.01
			(0.03)	(0.02)	(0.01)	(0.04)	(0.02)	(0.01)
	Medium	0.1	0.02	0.02	0.01	0.03	0.03	0.01
			(0.11)	(0.08)	(0.05)	(0.12)	(0.09)	(0.05)
		0.3	0.01	0.01	0.00	0.02	0.02	0.01
			(0.08)	(0.05)	(0.03)	(0.10)	(0.07)	(0.04)
		0.8	-0.01	0.00	0.00	0.01	0.01	0.00
			(0.05)	(0.04)	(0.02)	(0.08)	(0.06)	(0.04)
	High	0.1	0.02	0.01	0.01	0.00	0.00	0.00
			(0.09)	(0.07)	(0.04)	(0.09)	(0.07)	(0.04)
		0.3	0.01	0.00	0.01	0.00	0.00	0.00
			(0.07)	(0.05)	(0.03)	(0.08)	(0.06)	(0.03)
		0.8	0.01	0.00	0.00	-0.01	0.00	0.00
			(0.04)	(0.03)	(0.02)	(0.06)	(0.04)	(0.02)

Table 3 .

 3 Bias induced by a misspecification of population incidence on penetrance estimate (500 families, π =0.5, selection of families on at least two affected individuals, 200 replicates)

	Bias on penetrance estimate
	by age (years)

* When this factor is equal to 1, the true population incidence value is specified and the bias is due to limited sample size (see table

2)

Table 4 .

 4 Bias induced by a misspecification of the index case (500 families, π =0.5, 200 replicates) : minimal number of affected individuals required for selection of families; Youngest: analysis performed with the youngest affected specified as the index. Random: index chosen at random among affected individuals known as carriers. At least 1 A: no index specified but likelihood conditioned on at least one carriers among affected individual, tested or not; At least 1 AT: likelihood conditioned on at least one carrier among affected and tested individuals.

	True	Choice of index	Bias on penetrance estimate
	penetrance		MNA = 3	MNA = 4
	value		Age 50	Age 70 Age 50	Age 70
	Low	Youngest	0.00	0.02	0.00	0.02
		Random	0.01	0.01	0.01	0.01
		At least 1 A	0.02	0.06	0.02	0.06
		At least 1 AT	0.01	0.03	0.02	0.03
	Medium	Youngest	-0.04	0.04	-0.03	0.05
		Random	-0.02	0.01	-0.01	0.01
		At least 1 A	-0.02	0.06	-0.01	0.09
		At least 1 AT	-0.01	0.06	0.00	0.07
	High	Youngest	-0.06	0.00	-0.06	0.01
		Random	-0.01	-0.02	-0.01	-0.01
		At least 1 A	-0.03	-0.06	-0.02	-0.04
		At least 1 AT	0.01	0.01	0.01	0.01

MNA

Table 5 .

 5 Power to detect penetrance heterogeneity in case of mixture (50% each) of high and low risk families using a stratification on the number of families (500 families, π =0.5, MNA=2, 400 replicates)

			Estimated penetrance			Power for detecting
	By age	Assuming	After stratification of families according	heterogeneity
		homogeneity		to the number of affected individuals
			< 3	≥ 3	< 4	≥ 4
	40	0.17	0.12	0.14	0.18	0.19	0.08
	70	0.66	0.51	0.58	0.71	0.75	0.10

Table 6 .

 6 Comparison between single trait and multiple trait estimation for 4 cancer localizations by age 70 years according to MNA (500 families, π =0.5, 200 replicates)

		True
		penetrance	Single trait estimation	Multiple trait estimation
	Tumor
		value
		MNA	2	3	4	2	3	4
	CRC (men)	0.40	0.40	0.40	0.40	0.39	0.40	0.41
	CRC
		0.30	0.28	0.24	0.21	0.30	0.31	0.31
	(women)
	END	0.28	0.27	0.22	0.18	0.26	0.27	0.28
	OVA	0.07	0.07	0.06	0.05	0.08	0.08	0.08
	URE	0.02	0.02	0.015	0.014	0.04	0.04	0.04

CRC: colorectal cancer, END: endometrial cancer, OVA: ovarian cancer, URE: cancer of the urologic tract

Table 7 .

 7 Power for detecting a parent-of-origin effect according to number of families and genotyping rate (MNA=2, 200 replicates)

	Penetrance values
	Number of	Genotyping rate
	according to		Power
	families	(π)
	origin of mutation
		0.2	0.61
	100
		0.5	0.64
	P40=.25 & .40	0.2	0.97
	200
	P70=.30 & .60	0.5	0.98
		0.2	1.00
	500
		0.5	1.00
		0.2	0.01
	100
		0.5	0.01
	P40=.30 & .35	0.2	0.05
	200
	P70=.40 & .50	0.5	0.04
		0.2	0.17
	500

using our first version of the GRL [START_REF] Alarcon | Estimating cancer risk in HNPCC by the GRL method[END_REF] . Among studies that used a retrospective likelihood for estimating cancer risk in Lynch syndrome [START_REF] Alarcon | Estimating cancer risk in HNPCC by the GRL method[END_REF][START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF][START_REF] Jenkins | Cancer risks for mismatch repair gene mutation carriers: a population-based early onset case-family study[END_REF][START_REF] Stoffel | Calculation of risk of colorectal and endometrial cancer among patients with Lynch syndrome[END_REF] , note that only Quehenberger et al. [START_REF] Quehenberger | Risk of colorectal and endometrial cancer for carriers of mutations of the hMLH1 and hMSH2 gene: correction for ascertainment[END_REF] used a multiple trait approach. This may partly explain the variability of estimates among studies.

The test for a parent-of-origin effect was shown to be very efficient with samples of moderate size when the difference in penetrance was large. It would be interesting to apply the method for testing such an effect in hereditary forms of diseases, such as breast and ovarian cancers associated to mutations of BRCA1 or BRCA2 or Lynch syndrome associated with mutations of the MMR genes, as this hypothesis has not been tested yet.

The two extensions of the GRL have also been implemented in the proband's phenotype exclusion likelihood (PEL) that estimates penetrance in case of single ascertainment [START_REF] Alarcon | PEL: an unbiased method for estimating age-dependent genetic[END_REF] . Both methods are two different options of the GENERISK software which can be obtained from the authors (bernard.bonaiti@inserm.fr, BONADONA@lyon.fnclcc.fr, nadine.andrieu@curie.net).

CONFLICT OF INTEREST

The authors declare no conflict of interest disease risk from pedigree data unselected for family history. Genet Epidemiol 2009; 33: 379-385.

Appendix: example of misspecification of the index case

In 2005, a woman was affected by colorectal cancer at age 52 years and died two years later from ovarian cancer (case 1). Nobody suspected Lynch syndrome at that time although one of her uncle, currently aged 70 years, was affected at age 45 years by colorectal cancer (case 2). In 2009, the brother of the latter, aged 71 years, also developed colorectal cancer (case 3). The sister of case 1 suspected a hereditary syndrome and asked for genetic counselling. As case 1 was not available any more and case 3 could be a sporadic case because of late-onset diagnosis, the geneticist proposed genetic testing of case 2. Therefore, case 2 was designed as the index case although case 3 was obviously the incident case that would have been the 'natural' index.

Why should this choice induce a bias in the analysis? As the GRL is conditioned not only on the phenotypes of family members, but also on the genotype of the index case, the choice of case 2 as the index case 'cancels' his contribution to the likelihood and replaces it by the contribution of case 3. As case 3 was affected at an older age than case 2, this tends to overestimate age of onset and to decrease the penetrance estimate.


Alejandro Vazquez-Martin

Cristina Oliveras-Ferraros

Sonia Del Barco

Begoña Martin-Castillo

Ph.D Javier A Menendez
email: jmenendez@iconcologia.net

The anti-diabetic drug metformin suppresses self-renewal and proliferation of trastuzumab-resistant tumor-initiating breast cancer stem cells

Keywords: Metformin, Trastuzumab, HER2, Stem Cells, Tumor-initiating cells

à 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

For decades, breast cancer (BC) initiation and progression has been regarded as a multistep process that is reflected by the progressive genetic alterations that drive the transformation of normal human cells into highly malignant derivatives [START_REF] Jones | Cancer stem cells: are we missing the target?[END_REF]. In recent years, the experimental demonstration of tumor-initiating cells (popularly known as cancer stem cells) rather supports tumor hierarchy as a fundamental concept in the biology of BC disease [START_REF] Al-Hajj | Prospective identification of tumorigenic breast cancer cells[END_REF]. According to the cancer stem cell hypothesis, primary BC contain subpopulations of replenishing stem-like cells (i.e. tumor-initiating cells) with the ability of both self-renew and give rise to phenotypically diverse progeny [START_REF] Fillmore | Human breast cancer cell lines contain stem-like cells that self-renew, give rise to phenotypically diverse progeny and survive chemotherapy[END_REF]. In other words, solely the tumor-initiating cells have the ability to proliferate extensively, give rise to differentiated cells and form new tumors. When using cytotoxics or radiation, the bulk BC cell population decreases, but the proportion of tumor-initiating cells is higher than before treatment, thus indicating that conventional therapeutics efficiently kill differentiated BC cells whereas BC stem cells -by their nature-are intrinsically refractory thereby allowing BC regrowth [START_REF] Kakarala | Implications of the cancer stem-cell hypothesis for breast cancer prevention and therapy[END_REF][START_REF] Li | Intrinsic resistance of tumorigenic breast cancer cells to chemotherapy[END_REF]. The fact that tumorinitiating cells display intrinsic resistance to many conventional treatments might therefore explain the current limitations of these therapeutic approaches in curing breast BC.

Although many novel molecularly targeted agents (e.g. monoclonal antibodies, tyrosine kinase inhibitors [TKIs], etc.) have been also designed to target rapidly proliferating BC cancer cells -so many tumor-initiating cells might also be relatively insensitive to these agents-we are beginning to accumulate evidence that anti-HER2 therapeutics can effectively target tumor-initiating cells derived from HER2overexpressing BC [START_REF] Li | Intrinsic resistance of tumorigenic breast cancer cells to chemotherapy[END_REF][START_REF] Magnifico | Tumor-initiating cells of HER2-positive carcinoma cell lines express the highest oncoprotein levels and are sensitive to trastuzumab[END_REF]. In analogy with genes such as Notch, Hedgehog, Wnt and ras, initially identified as oncogenes but shown subsequently to be actively involved in selfrenewal and survival of tumor-initiating cells, HER2 may also constitute a pivotal molecule in BC stem cells. The notion that HER2 signaling can drive tumor-initiating cells in certain HER2-positive BC is supported by the following findings: 1.) HER2 overexpression in multiple BC cell lines significantly increases the fraction of cells positive for ALDH1 (aldehyde dehydrogenase) [START_REF] Korkaya | HER2 regulates the mammary stem/progenitor cell population driving tumorigenesis and invasion[END_REF] -a marker of normal and cancerous human mammary epithelial cells with stem/progenitor properties-and enhances the tumor-initiating component that drives tumorigenesis, invasion and metastasis [START_REF] Ginestier | Stemming resistance to HER-2 targeted therapy[END_REF] To define the sensitivity of the sub-population of tumor-initiating (BC stem) cells to current and novel anti-HER2 BC therapies; and 2.) To discover agents and strategies that can target HER2-positive BC and HER2-positive BC relapse at their apparent source of tumor-initiating (BC stem) cells.

Culturing BC cells in so-called "mammosphere culture" has been proposed as a reliable method for enriching sub-populations of putative BC tumor-initiating stem cells [START_REF] Grimshaw | Mammosphere culture of metastatic breast cancer cells enriches for tumorigenic breast cancer cells[END_REF]. Mammosphere culture could provide a highly appropriate model for studying the sensitivity of the tumor-initiating stem cells within a heterogeneous BC cell population to therapeutic agents. Thus, a comparison of efficacy using mammospheres from BC cells lines or mammospheres cultured from individual pleural effusions as targets would allow evaluation of the consistency of effectiveness of a specific therapy.

Taking advantage of the ability of BC stem cells to form multicellular "microtumors" in non-adherent and non-differentiating conditions (i.e. "mammospheres" [START_REF] Dontu | In vitro propagation and transcriptional profiling of human mammary stem/progenitor cells[END_REF][START_REF] Dontu | Role of Notch signaling in cell-fate determination of human mammary stem/progenitor cells[END_REF]) and considering novel data showing the unexpected ability of the anti-diabetic drug metformin to preferentially kill BC stem cells [START_REF] Hirsch | Metformin selectively targets cancer stem cells, and acts together with chemotherapy to block tumor growth and prolong remission[END_REF], we recently decided to explore if metformin may alter self-renewal and proliferation of tumor-initiating stem cells in BC in vitro models representative of Tzb-responsive and Tzb-refractory HER2-gene amplified BC disease. We have assessed also, at the level of mammosphere-initiating cells, the relevance of well-recognized BC stem cell markers (i.e. CD44) in the phenomena of both acquired and de novo refractoriness to the anti-HER2 monoclonal antibody Tzb. Finally, we have explored the nature of the (anti-BC stem cell) interaction (i.e. antagonism, addition, and synergism) between metformin and Tzb in Tzb-naïve and Tzb-refractory HER2-overexpressing BC models.

MATERIALS and METHODS

Materials. Trastuzumab (Herceptin™) was provided by Pharmacy Department, Hospital Universitario de Girona Dr. Josep Trueta (Girona, Spain), solubilized in bacteriostatic water for injection containing 1.1% benzyl alcohol (stock solution at 21 mg/ml), stored at 4°C and used within one month. Metformin hydrochloride (1,1dimethylbiguanide hydrochloride) was purchased from Sigma Chemical (St. Louis, MO, USA), solubilized in phosphate-buffered saline (PBS; stock solution at 1 mol/L), stored at -20 o C and used within two weeks.

Culture Conditions. SKBR3, MDA-MB-468 and MDA-MB-231 human breast cancer cells were obtained from the American Type Culture Collection (ATCC) and they where routinely grown in Improved MEM (IMEM; BioSource International; Invitrogen S.A., Barcelona, Spain) supplemented with 5% fetal bovine serum (FBS) and 2 mM L-Glutamine. SKBR3 TzbR cells (POOL1), an in vitro model with acquired autoresistance to Tzb that grows in the continuous presence of 100 g/mL Tzb, has been developed in our laboratory as previously described [START_REF] Vazquez-Martin | Autophagy facilitates the development of breast cancer resistance to the anti-HER2 monoclonal antibody trastuzumab[END_REF]. JIMT-1 human breast cancer cell line was established at Tampere University from a ductal carcinoma pleural metastasis of a 62-year-old-BC patient who did not respond to Tzb [START_REF] Tanner | Characterization of a novel cell line established from a patient with Herceptin-resistant breast cancer[END_REF]. JIMT-1 cell line is available from the German Collection of Microorganisms and Cell Cultures (http://www.dsmz.de/). JIMT-1 cells were grown in F-12/DMEM (1:1) supplemented with 10% FBS and 2 mM L-glutamine. Cells were maintained at 37 o C in a humidified atmosphere of 95% air and 5% CO 2 . Cells were screened periodically for Mycoplasma contamination.

Mammosphere culture. Mammospheres were generated from single cells of SKBR3, SKBR3 TzbR, and JIMT-1 cell lines seeded at 10 3 cells/cm 2 in six-well ultralow attachment plates (Corning Inc.). Sphere medium (F-12/DMEM containing 5 mg/mL insulin, 0.5 mg/mL hydrocortisone, 2% B27 [Invitrogen Ltd.], and 20 ng/mL epidermal growth factor) was prepared as originally described by Dontu et al. [START_REF] Dontu | In vitro propagation and transcriptional profiling of human mammary stem/progenitor cells[END_REF].

Mammosphere-forming efficiency (MSFE). MSFE was calculated as the number

of spheres-like structures (large diameter > 50 m) formed in 7 days divided by the original number of cells seeded and expressed as percentage means ( SD) Flow cytometry. C ell cultures growing in regular medium (or supplemented with Tzb and/or metformin as specified), were washed once with phosphate-buffered saline (PBS) and then harvested with 0.05% trypsin/0.025% EDTA into single cell suspensions. Detached cells were washed with PBS containing 1% FBS and 1% penicillin/streptomycin (wash buffer), counted, and resuspended in the wash buffer (10 6 cells/100 l). Combinations of fluorochrome-conjugated monoclonal antibodies obtained from BD Biosciences (San Diego, CA, USA) against human CD44 (FITC; cat.#555478) and CD24 (PE; cat.#555428) or their respective isotype controls were added to the cell suspension at concentrations recommended by the manufacturer and incubated at 4 o C in the dark for 30-40 min. Labeled cells were washed in the wash buffer to eliminate unbound antibody, then fixed in PBS containing 1% paraformaldehyde, and then analyzed no longer than 1 h post-staining on a BD FacScalibur (BD Biosciences).

Statistical analysis.

Results are presented as the mean SD for at least three repeated individual experiments for each group. Analysis was done using XLSTAT 2010 (Addinsoft™). P 0.01 was considered statistically significant.

RESULTS

Metformin is as efficient as Tzb at reducing MSFE in Tzb-naïve SKBR3 BC cells. The mammosphere culture technique was employed to assess whether Tzb and metformin, as single agents, may significantly reduce MSFE in Tzb-naïve SKBR3 BC cells (Fig. 1A). A clinically relevant concentration of Tzb was sufficient to significantly decrease MSFE in HER2-overexpressing SKBR3 cells ( 35% MSFE reduction at 10 g/mL Tzb; Fig. 1B). Similarly, a 35% MSFE reduction was obtained when sphere cultures of SKBR3 cells were carried out in the presence of 100 mol/L metformin (Fig. 1B). Indeed, a dose-dependent reduction of MSFE was clearly observed when SKBR3 cells were cultured in sphere medium supplemented with graded concentrations of metformin (up to 75% reduction at 1000 mol/L metformin; Fig. 1A,B). Remarkably, metformin treatment not only reduced MSFE but caused further morphological changes to SKBR3 mammospheres. Thus, SKBR3 mammospheres formed in the presence of 50 and 100 g/mL metformin exhibited a significant disintegration of the spherical mammospheres structures observed in untreated SKBR3 control cells (Fig. 1A). When sphere cultures of SKBR3 cells were carried out in the presence of concurrent Tzb and metformin, we observed mostly additive interactions of both agents when analyzing fold-reductions of MSFE in SKBR3 cells (Fig. 1C). remained largely unaltered when SKBR3 TzbR sphere cultures were performed in the presence of Tzb doses as high as 100 g/mL (Fig. 2B,C).

Acquisition of auto-resistance to

Metformin interacts synergistically with Tzb at reducing MSFE in HER2-

overexpressing BC cells with acquired resistance to Tzb. In light of these findings, we sought to assess whether metformin could reverse acquired autoresistance to Tzb by targeting "Tzb-resistant" HER2-positive BC stem/progenitor cells.

Of note, Tzb-refractory mammosphere-initiating cells were exquisitely sensitive to metformin. 100 mol/L metformin reduced MSFE by about 50% while 1000 mol/L metformin completely suppressed mammosphere formation in SKBR3 TzbR cells (Fig. 2A,B). Moreover, synergistic interactions between Tzb and metformin were found throughout the entire range of metformin concentrations (Fig. 2C). Thus, the ability of 500 mol/L metformin to decrease SKBR3 TzbR MSFE increased from about 5 times when used as single-agent to more than 17 times when combined with Tzb. JIMT-1 cells cultured in sphere medium for 7 days formed irregular and somewhat loosely-adhered clumps of cells (Fig. 3A). Although treatment with single-agent Tzb reduced both the number and the size of JIMT-1 mammosphere-like structures, Tzbinduced reduction of JIMT1 MSFE was not statistically significant (Fig. 3B,C). These findings support further the notion that refractoriness to Tzb in HER2-positive BC cells can be explained in terms of BC stem/progenitor populations intrinsically resistant to Tzb. In this scenario, treatment with single-agent metformin was sufficient to drastically reduce MSFE in Tzb-refractory JIMT-1 cells ( 65% MSFE reduction at 100 mol/L metformin; Fig. 3B,C). Perhaps more importantly, concurrent combination of Tzb and metformin worked synergistically to totally suppress the ability of JIMT-1 cells to form mammospheres. Thus, the ability of 1000 mol/L metformin to decrease MSFE increased from about 5 times when used as single-agent to more than 22 times when combined with Tzb (Fig. 3C).

Metformin interacts synergistically with

Metformin interacts synergistically with Tzb to suppress the CD44 pos CD24 neg/low stem/progenitor immunophenotype in Tzb-refractory JIMT-1 BC cells. Given that Tzb-refractory HER2-overexpressing JIMT-1 cells constitutively express significant amounts of the transmembrane receptor CD44 [START_REF] Pályi-Krekk | Hyaluronan-induced masking of ErbB2 and CD44-enhanced trastuzumab internalisation in trastuzumab resistant breast cancer[END_REF] -a cardinal feature of BC stem cells-we sought to examine: 1.) Whether unsorted populations of JIMT-1 differentiated cultures were naturally enriched with the stem/progenitor CD44 pos /CD24 neg/low progeny; and 2.) Whether concomitant exposure to Tzb and metformin significantly affected the BC stem CD44 pos /CD24 neg/low immunophenotype. When using JIMT-1 cells from cultures passages #20 to #30, flow cytometry analyses clearly revealed a significant subpopulation of about 70% of JIMT-1 cells naturally exhibiting the CD44 pos /CD24 neg/low BC stem profile (Fig. 3D). Rather than exhibiting a luminal-type classification, the unique distribution of CD24 & CD44 markers in HER2-gene amplified JIMT-1 cells somewhat resembled that observed in CD44 pos /CD24 neg/low -enriched MDA-MB-231 mesenchymal BC cells (Fig. 3A). Then, we assessed distribution of the cell surface markers CD24 and CD44 in JIMT-1 cells following treatment (48 h) with two concentrations of metformin (i.e. 100 and 1000 mol/L) in the absence or presence of an optimal dose (100 g/mL) of Tzb (Fig. 3D).

Whereas single-agent Tzb failed to significantly reduce the number of JIMT-1 cells exhibiting the putative stem-like CD44 pos CD24 neg/low immunophenotype, CD44 pos CD24 neg/low cell numbers were reduced by ~30% upon treatment with 1000 g/mL metformin. Concurrent treatment with Tzb and 1000 mol/L metformin reduced further, in a synergistic fashion, CD44 pos CD24 neg/low cell numbers (from 68% in untreated control cells to 37% in Tzb plus metformin co-treated cells) to concurrently increase the CD44 pos CD24 pos fraction (from 32% in untreated control cells to 62% in Tzb plus metformin co-treated cells) (Fig. 3D).

To exclude that synergistic reduction of MSFE efficiencies in Tzb-refractory SKBR3

TzbR and JIMT-1 cells would simply reflect a generalized reduction in cellular viability, we replicated all the experimental groups in adherent conditions and cell growth was monitored by MTT assay. At the doses employed in sphere cultures, Tzb and metformin were ineffective at inhibiting cell growth of SKBR3 TzbR and JIMT-1 differentiated cultures routinely grown in plastic (data not shown).

DISCUSSION

Hirsch et al. have recently demonstrated that tumor-forming, self-renewing BC stem cells appear to exhibit an exacerbated sensitivity to the anti-diabetic drug metformin [START_REF] Hirsch | Metformin selectively targets cancer stem cells, and acts together with chemotherapy to block tumor growth and prolong remission[END_REF]. In their hands, concurrent treatment with metformin (100 or 300 mol/L) and the cytotoxic doxorubicin was found to reduce tumor mass and prevent relapse much more effectively than either drug alone in a xenograft mouse model. Here, using three with graded doses of metformin, it was relevant also that metformin caused gross morphological changes in mammosphere structure. It is reasonable to suggest that, during formation of multicellular "microtumors", metformin treatment may actively interfere with the basic processes of differentiation and/or apoptosis in the stem/progenitor cell population of Tzb-naïve, HER2-overexpressing BC cells.

HER2
Second, our current findings provide a novel scenario for the molecular understanding of Tzb resistance in HER2 gene-amplified BC (Fig. 4). The acquisition of Tzb resistance upon chronic exposure to Tzb can be explained in terms of emerging "Tzb-resistant" CD44-positive BC stem/progenitor cells because, in our hands, SKBR3

TzbR mammosphere-initiating cells did not respond to Tzb exposure. Our current findings strongly suggest that during acquisition of auto-resistance to Tzb (i.e.

secondary resistance to Tzb), HER2-overexpressing luminal BC cells appear to undergo a transition to basal-like phenotypes characterized by a significant enrichment in cell populations overexpressing the stem-cell marker CD44. This CD44-related scenario for Tzb refractoriness was even more apparent in HER2-gene amplified BC cells with de novo resistance to multiple HER-targeted therapies. JIMT-1 BC cells not only exhibit primary resistance to Tzb but demonstrate further cross-resistance to several HER2inhibiting drugs including the HER-dimerization inhibitor antibody pertuzumab (2C4)

and the small molecule tyrosine kinase inhibitors (TKIs) CI1033 (canertinib) and ZD1839 (gefitinib). Preliminary experiments in our laboratory demonstrated that JIMT-1 cells display also cross-resistance to the dual HER1/HER2 TKI lapatinib (data not shown). Given that Tzb treatment has been shown to fail at decreasing the fraction of JIMT-1 cells positive for ALDH1 [START_REF] Korkaya | HER2 regulates the mammary stem/progenitor cell population driving tumorigenesis and invasion[END_REF], it has been suggested that BC stem cell-like phenotypes may underlie Tzb resistance in JIMT-1 cells [9]. In this regard, we have validated the notion that Tzb-refractory HER2-positive JIMT-1 cell cultures contain a significant population bearing the CD44 pos CD24 neg/low phenotype initially described by Al-Hajj et al. as a feature of BC stem cells that associates with the expression of basal/mesenchymal or myoepithelial markers [START_REF] Al-Hajj | Prospective identification of tumorigenic breast cancer cells[END_REF]. Because earlier studies have repeatedly confirmed that all BC cell lines with significant CD44 pos /CD24 neg/low subpopulation belong to the basal type of BC and display the phenotype of cancer cells that have undergone epithelial to mesenchymal transition (EMT) [START_REF] Sheridan | CD44+/CD24-breast cancer cells exhibit enhanced invasive properties: an early step necessary for metastasis[END_REF][START_REF] Honeth | The CD44+/CD24-phenotype is enriched in basal-like breast tumors[END_REF], the unique distribution of CD24 and CD44 cell surface markers in HER2-overexpressing JIMT-1 cells may suggest that the induction of EMT in HER2-positive BC cells may generate cells with stem-like features, such as the CD44 pos /CD24 neg/low phenotype and the ability to grow in the presence of HER2 inhibitors. In this regard, it is noteworthy that JIMT-1 cells failed to produce discrete and correctly organized spherical structures but rather formed somewhat loosely adhered clumps of cells when cultured in sphere medium for 7 days, a mammosphere-like architecture that has been observed in sphere cultures of CD44 pos /CD24 neg/low -enriched HER2-negative MDA-MB-231 mesenchymal BC cells [START_REF] Grimshaw | Mammosphere culture of metastatic breast cancer cells enriches for tumorigenic breast cancer cells[END_REF].

Third, mammosphere-initiating cells from both CD44-overexpressing SKBR3 TzbR cells with acquired auto-resistance to Tzb and CD44-overexpressing JIMT-1 cells with de novo resistance to Tzb retained an exquisite and somewhat exacerbated sensitivity to the inhibitory effects of metformin. It has previously been shown that mammosphere number reflects stem cell self-renewal, whereas mammosphere size reflects progenitor cell proliferation [START_REF] Dontu | In vitro propagation and transcriptional profiling of human mammary stem/progenitor cells[END_REF][START_REF] Dontu | Role of Notch signaling in cell-fate determination of human mammary stem/progenitor cells[END_REF]. Our current findings reveal that Tzb treatment may affect stem/progenitor cell proliferation in Tzb-refractory HER2-gene amplified BC because it significantly reduced mammosphere size. The incapability of Tzb to significantly decrease mammosphere cell number strongly suggests that Tzb efficacy in HER2overexpressors rather relates to its ability to block HER2-driven stem/progenitor cell self-renewal. Of note, metformin co-treatment synergistically worked with Tzb to suppress self-renewal and proliferation of Tzb-refractory tumor-initiating cells In this scenario, it would be tempting to suggest that metformin-targeted crucial survival mechanisms may underlie also Tzb resistance in tumor-initiating cells because concurrent combination of metformin with Tzb synergistically reduced MSFE as well as mammosphere size in Tzb-refractory SKBR3 TzbR and JIMT-1 cells. Given that the percentage of CD44 pos CD24 neg/low tumor-initiating cells has been shown to significantly correlate with MSFE in BC [START_REF] Li | Intrinsic resistance of tumorigenic breast cancer cells to chemotherapy[END_REF] and that co-treatment with metformin and Tzb synergistically reduced the percentage of JIMT-1 cells exhibiting the CD44 pos CD24 neg/low stem/progenitor immunophenotype, it would be relevant to elucidate whether, in the presence of metformin, "Tzb-resistant" BC stem/progenitor cells may regain sensitivity to Tzb.

In summary, we present the first evidence that the anti-diabetic drug metformin synergistically interacts with molecularly targeted agents such as Tzb to suppress selfrenewal (measured as MSFE) and proliferation (evaluated as mammosphere size) of tumor-initiating (BC stem) cells of HER2-positive carcinoma cells showing acquired and de novo resistance to Tzb. The therapeutic plasma levels of metformin in clinical medicine are 0.465 to 2.5 mg/mL, whereas we observed statistically significant reductions in MSFE and mammosphere size in 50 to 1000 g/mL metformin, a minimum 3-fold excess over the recommended therapeutic levels. Although the range of concentrations employed in our current approach were significantly lower than those used in earlier studies exposing cultures of differentiated cancer cells to metformin concentrations ranging from 1 to 100 mmol/L, forthcoming studies should validate whether concentrations of metformin capable of efficiently eliminating mammosphereinitiating cells might be also attained during breast cancer treatment. Although it has been reported that metformin can accumulate in tissues at concentrations several-fold higher than those in blood [START_REF] Wilcock | Accumulation of metformin by tissues of the normal and diabetic mouse[END_REF], caution must be therefore applied when extrapolating our current in vitro results into clinical practice. Nevertheless, it is reasonable to suggest that, in combination with conventional chemotherapy and HER2-targeting

  ; 2.) Pioneering findings by Li et al demonstrated that, conversely to chemotherapy treatment, the epidermal growth factor receptor (EGFR)/HER2 dual TKI lapatinib, decreases the subpopulation of chemotherapy-resistant BC initiating cells (measured by calculating the percentage of cells with the CD44 pos /CD24 neg/low stem/progenitor immunophenotype and the ability to form mammospheres in vivo as an indication of self-renewal) [5], and; 3.) Tumor-initiating cells of HER2-positive BC have been shown to express the highest HER2 levels and been sensitive to the anti-HER2 monoclonal antibody trastuzumab [Tzb; Herceptin ] [6]. Intriguingly, recent observations have suggested also that clinical resistance to Tzb and other HER2-targeted therapies may be driven by stem-like BC cells within HER2-positive BC. Indeed, a variety of possible mechanisms of escape from Tzb involves many of the same markers that have been implicated in the biology of cancer stem cells (e.g. loss or blockade of the Tzb binding site at the extracellular domain of HER2 upon overexpression of the stem cell marker CD44, signaling through alternative receptor pathways via the stem cell markers CXCR4 and 1 integrin, activation of Notch signaling to promote Tzb resistance via enhanced self-renewal/proliferation of progenitor cells or induction of epithelia-tomesenchymal transition, etc.) [Reviewed in 9]. Therefore, a major challenge now is: 1.)

  Tzb can be explained in terms of "Tzbresistant" HER2-positive BC stem/progenitor cells. Given the remarkable heterogeneity of BC disease, it has been suggested that the characteristics of HER2positive BC tumor-initiating cells may evolve under the selection pressure of Tzb to drive the emergence of "Tzb-resistant" HER2-positive BC stem/progenitor cells [9]. To evaluate this hypothesis, we took advantage of SKBR3 TzbR cells (POOL1), which were obtained in our laboratory by continuously culturing HER2-dependent SKBR3 parental cells in the presence of high-doses of Tzb for more than 10 months[START_REF] Vazquez-Martin | Autophagy facilitates the development of breast cancer resistance to the anti-HER2 monoclonal antibody trastuzumab[END_REF]. This Tzbconditioned SKBR3 derivative optimally grows in the presence of 100 g/mL Tzb.Unlike Tzb-sensitive SKBR3 parental cells, which formed discrete and spherical mammospheres, SKBR3 TzbR cells cultured in sphere medium formed larger & irregular mammosphere-like structures (Fig.2A). To evaluate whether changes in gross mammosphere structure related to changes in the expression of the stem cell markers involved also in Tzb resistance, we employed flow cytometry to examine cell surface distribution of CD24 and CD44 markers before and after acquisition of Tzb autoresistance. Tzb-naïve SKBR3 parental cells likewise exhibited a significantly enriched CD24 pos /CD44 neg ( 90%) cell population, which is consistent with the luminal-type classification of SKBR3 BC cells (Fig.2A, left). Interestingly, 55% of the SKBR3 TzbR cell population displayed a CD24 pos /CD24 pos immunophenotype, which to some extent resembled the CD24/CD44 distribution that can be observed in highly-metastatic MDA-MB-468 basal-like BC cells (Fig.2A, right). Moreover, mammosphere-initiating cells within SKBR3 TzbR cell populations did exhibit resistance to Tzb as MSFE

  Tzb to reduce MSFE in JIMT-1 BC cells, an in vitro model of refractoriness to HER2-targeted therapies ab initio. To validate further the notion that metformin may attack Tzb-resistant tumorinitiating stem cells, we took advantage of the de novo Tzb resistance in a cell line (JIMT-1) isolated from the pleural fluid of a HER2-positive BC patient with progressive disease on Tzb [15]. Unlike Tzb-sensitive HER2-overexpressing SKBR3 cells and similarly to SKBR3 TzbR cells with acquired auto-resistance to Tzb, Tzb-refractory

  -gene amplified BC in vitro models for Tzb-sensitivity (SKBR3), acquired Tzb auto-resistance (SKBR3 TzbR) and de novo Tzb refractoriness (JIMT-1), and employing mammosphere-forming efficiency & mammosphere size as indicators of stem selfrenewal and progenitor cell proliferation, respectively, we provide for the first time preclinical evidence for the efficacy of the anti-diabetic drug metformin (50-100 mol/L) to efficiently target Tzb-resistant tumor-initiating stem cells in HER2 overexpressors. We also confirm recent findings by Magnifico et al. showing that tumor-initiating cells of HER2-positive carcinoma cell lines are sensitive to the anti-HER2 monoclonal antibody Tzb [6]. Therefore, our findings corroborate and expand further Hirsch and Magnifico pioneer's findings to the combination of metformin with Tzb against HER2positive tumor-initiating BC stem cells. First, we confirmed the distinctive ability of the anti-diabetic drug metformin to attack BC tumor-initiating stem cells within a heterogeneous population of BC cells overexpressing HER2. It should be noted that metformin concentrations as low as 50 mol/L were sufficient to mimic the therapeutic effects of Tzb against mammosphereinitiating cells in Tzb-naïve SKBR3 cells. Although the remarkable potency to eliminate stem/progenitor cells in sphere cultures of HER2 overexpressors became obvious when analyzing dose-dependent decreases in MSFE of SKBR3 sphere cultures supplemented

Figure 1 .

 1 Figure 1. Metformin is as efficient as Tzb at inhibiting mammosphere

Figure 2 .

 2 Figure 2. Metformin and Tzb interact synergistically to suppress self-

Figure 4 .

 4 Figure 4. Metformin and Tzb-resistant tumor-initiating BC stem cells: A

agents, metformin co-treatment may add an important therapeutic benefit in debulking and in targeting tumor-initiating cells of HER2-overexpressing BC [20,21]. Although it remains to be elucidated which metformin-sensitive molecular mechanisms control the maintenance and survival of Tzb-resistance HER2-gene amplified tumor initiating cells, crucial evidence of the clinical efficacy of MET could be obtained from "proof-ofprinciple" studies with neoadjuvant regimens including neoadjuvant chemotherapy (to target bulk BC cells) and Tzb with MET (both targeting self-renewal and proliferation of tumor-initiating cells) in women diagnosed with HER2-positive primary BC [22].

Acknowledgments

Alejandro Vazquez-Martin is the recipient of a "Sara Borrell" post-doctoral contract (CD08/00283, Ministerio de Sanidad y Consumo, Fondo de Investigación Sanitaria -FIS-, Spain). Work at Menendez' Laboratory is supported by Instituto de Salud Carlos III (Ministerio de Sanidad y Consumo, Fondo de Investigación Sanitaria -FIS-, Spain, Grants CP05-00090 and PI06-0778 and RD06-0020-0028), the Fundación Científica de la Asociación Española Contra el Cáncer (AECC, Spain), and by the Ministerio de Ciencia e Innovación (SAF2009-11579, Plan Nacional de I+D+ I, MICINN, Spain).

Breast Cancer Research and Treatment

The anti-diabetic drug metformin suppresses Self-renewal and proliferation of trastuzumab-resistant Tumor-initiating breast cancer stem cells


Pablo Cortes

Pablo Cortés

Luis Onieva

José Guadix

Optimizing and simulating the assembly line balancing problem in a motorcycle manufacturing company

Keywords: ASSEMBLY LINE BALANCING, LEAD-TIME REDUCTION, SIMULATION motorcycle, assembly line, lead time, simulation

Optimizing and simulating the assembly line balancing problem in a motorcycle manufacturing company: a case study

Introduction

Automobile manufacturing is the industry sector that makes greatest use of assembly lines. The problem of design and balancing of assembly lines is difficult to be solved due to its combinatorial nature [START_REF] Wee | Assembly line balancing as generalized bin packing[END_REF] and also due to the large amount of constraints and tasks to be carried out in real life problems.

An assembly line is a manufacturing process in which parts are added to a product in a sequential manner using optimally planned logistics to create a finished product much faster than with handcrafting-type methods. Essentially, the main aspects of assembly lines are the number of stations and the production rate. In the case we are dealing with in our paper, every worker is supposed to have the same amount of time for carrying out the operations that are associated to his/her post, C.

Consequently, 1/C is the production rate better performance for the whole of experiments, which should be sufficient to meet the required demand. Operations (or tasks) are assigned to the stations where the workers assemble the parts and subcomponents, depending on the precedence constraints between operations, the fixed stations constraints where an operation may have to be carried out, and the incompatible tasks constraints for operations that cannot be carried out in the same station. Assembly lines can also include many other specific constraints such as linked tasks, cumulated restrictions, excluded stations, station type, and minimal/maximum distance, amongst others. An interesting web page addressing most of the well-known problems in assembly lines can be accessed in the site www.assembly-line-balancing.de. This paper discusses manufacturing assembly line performance improvement of a company that produces motorcycles named MH, whose main factory is located in Spain. It is the most southern automobile industry company in Europe and the majority of its production (more than 75%) is destined for export. Currently, the most outstanding product of the company is the FURIA XP6 model (Figure 1), which is the model we analyse in this paper.

FIGURE 1

This paper focuses on the assignment of operations to the assembly line, comparing the previous scheduling in the company with diverse proposals. We address a method based on the traditional rank positional weight heuristic, which is adapted to the specific conditions of motorcycle assembly when two people must be always working on the vehicle at each station, and a neighbourhood search heuristic, which provided the better quality results. This neighbourhood search heuristic In a general scope the assembly line problem has been widely studied in scientific literature, being surveyed by several authors [START_REF] Rekiek | State of art of optimization methods for assembly line design[END_REF][START_REF] Becker | A survey on problems and methods in generalized assembly line balancing[END_REF][START_REF] Scholl | State-of-the-art exact and heuristic solution procedures for simple assembly line balancing[END_REF][START_REF] Boysen | A classification of assembly line balancing problems[END_REF] as recent examples). Focusing on real case studies, an interesting paper [START_REF] Lapierre | Balancing assembly lines: an industrial case study[END_REF]) addresses a case study in the Canadian industry of appliances. The paper considers a two-sided conveyor with different speeds and grouping characteristics that make it different from other traditional cases. The problem was solved by implementing a priority-based heuristic developed in MsAccess. However, there are not many research papers dealing with the specific problem of motorcycle production. [START_REF] Pastor | Assembly line balancing with incompatibilities and bounded workstation loads[END_REF] studied the assembly line balancing problem with incompatibilities and bounded workstation loads in a motorcycle factory. They presented the case of a single-sided line with a given number of stations and different types of assignment restrictions: some tasks are required to be assigned to the same station and, thus, they may be combined to a single task; some tasks can only be performed at the left hand or the right hand side of a large, irremovable workpiece; and some tasks need to be executed at the top or at the bottom of the workpiece. The specific characteristics of the current assembly line imply a number of specific constraints that must be taken into account when considering their paper, differing from the case study we are addressing here. More recently, some of the same authors [START_REF] Corominas | Balancing assembly line with skilled and unskilled workers[END_REF] have presented the process of rebalancing in a motorcycle assembly line with the aim of reducing the number of temporary workers. This problem considers workers with different skills that must be correctly assigned to stations in order to achieve the maximum productivity. Their paper differs in this respect from ours where all the workers are stable in the company staff. Additionally, we can assume that resources are homogeneous, and hence the task duration is not depending on its workstation allocation. Some of the same authors have continued researching on motorcycle assembly lines cases; it is the case of the alternative subgraphs assembly line balancing problem recently addressed in [START_REF] Capacho | ASALBP: the alternative subgraphs assembly line balancing problem[END_REF], or the papers by [START_REF] Andrés | Balancing and scheduling tasks in assembly lines with sequence-dependent setup times[END_REF], and [START_REF] Pastor | An improved mathematical program to solve the simple assembly line balancing problem[END_REF]. However, our paper relates more to [START_REF] Akagi | A method for assembly line balancing with more than one worker in each station[END_REF] and [START_REF] Lutz | Development of operator assignment schedules: a DSS approach[END_REF] where the performance of workers on a specific task is the same. The difference from our work relies on the specific characteristic of a motorcycle assembly line where exactly two persons have to work at the same time at each station.

The rest of the paper follows with the second section, which includes the description of the assembly line and the operations that are needed to produce a FURIA XP6 motorcycle model. The third section details the previous scheduling in the assembly line in the company, as well as the proposals that were considered.

Expected results from a theoretical point of view are also shown. The fourth section is dedicated to the simulation model and the results that the simulation model provided using ARENA software. Finally, the main aspects are reviewed in the conclusion.

Description of the assembly line for FURIA XP6 model

The main parts of a motorcycle are: the chassis, which constitutes the backbone of the motorcycle and is where the rest of the components are assembled; the clutch-release yoke, which joins the handlebars to the front wheel and acts as front suspension; a tipping device to join the rear wheel and the rear muffler; the petrol deposit; the engine, which is joined to the dumper and to the rear wheel by a chain or belt with the The morphology of motorcycles allows two operators to work on a vehicle at the same time but on separate tasks, or while collaborating on a more complex task. 57 operations are carried out for the FURIA XP6 model and are described in Table 1, where the average time and standard deviation have been measured for each task. In our case, we need to consider precedence constraints between operations, some fixed stations constraints where an operation has to be carried out (e.g. the chassis assembly), and incompatible tasks constraints for operations that cannot be carried out at the same station. The table also indicates when an operation must be carried out by both operators (B), and when by the right (R) or left (L) one. The table states the operations as they were carried out in the initial scenario.

TABLE 1

We have to note related to Table 1, that we carried out a checking analysis trying to identify which tasks were forced to be undertaken by "both" operators with the objective to relax this restriction in an hypothetical re-assembly. We proposed to assign operation 7 (firstly assigned to B) to the left operator due to assembly reasons, and tasks 14, 22, 48, and 51 that were initially conditioned to both operators to one of them (indistinctively to the left or right operator). These proposals were checked during the Christmas planned maintenance stop and then approved by the General Manager after checking its feasibility.

Once we have calculated the average times of each operation, t i , for N operations, we can calculate the total assembly time, T, and given a number of stations, K, and a lead time, C, we can calculate the inactivity time

T C K I - ⋅ = .
Thus, the objective of our problem is to minimize the lead time, C.

Figure 3 depicts workers on the assembly line carrying out operations detailed in Table 1.

... , { 2 1 N J J J
the set operations to be carried out on the assembly line. The execution time for the operations is known from Table 1, and is given by

i t , N i ,...... 1 =
. The precedence constraints among operations are also known. As previously stated, the unique feature of motorcycle assembly lines is that two operators can carry out operations at the same time, or cooperate to carry out a more complex operation. The goal is to assign operations to stations with the objective of minimizing the lead time.

Previous situation

The initial scheduling of the assembly line is described in Table 2. Table 2 shows the operations that have to be carried out by two operators simultaneously, due to the special requirements of the task. The other operations can be interchangeable between the right and left operator.

The table also includes waiting times. Idle times can result due to: (i) worker

A is waiting for worker B to make a joint operation or (ii) worker A is waiting for worker B because worker A's operation relies on worker B's operation being completed.

Modified Rank Positional Weight Heuristic (MRPWH)

First, we tried a new scheduling given by a heuristic based on the Rank Positional Weight Heuristic (RPWH) adapted to the unique situation of two operators working on the same vehicle on an assembly line. RPWH is an old reputed heuristic method [START_REF] Helgeson | Assembly Line Balancing Using Rank Positional Weight Technique[END_REF]. We detail our modification of RPWH below in the following rules:

• A weight is assigned to each operation i J according to

∑ ∈ + = D i i i i t t p
. The weight is calculated as the operation i J time plus the time of all the operations that must follow operation i J , that are stated as set D in the formula. • Operations are ordered in decreasing order of weight, safeguarding the precedence relationship. This configures an assignment list. • An approximate lead time is estimated as the sum of the operations' average time, divided by the number of stations and then multiplied by two due to the number of operators (in our case: 8). • Operations are assigned to operators following the order of the assignment list and according to the position requirements (left, right or both) up to satisfy the lead time (the nearest possible value to the lead time). First, in case of position requirement, the constraint is satisfied by assigning the operation to the specified operator. Second the rest of operations are assigned to left or right operator attending to his/her work load in the station and with the objective of minimizing the idle time.

• The station's total time is calculated for the four stations.

• The lead time is calculated as the maximum of the four stations.

Table 3 shows the scheduling proposal according to the MRPWH technique. As can be seen, the lead time is reduced from 0:12:49 to 0:07:29, which is approximately 41% of the time saved, after applying MRPWH.

Multi-Started Neighbourhood Search Heuristic (MSNSH)

We developed a novel multi-started neighbourhood search heuristic, with the aim of providing a more balanced assembly line.

FIGURE 4

Given a feasible solution, x i , the quality of the feasible solution, ( )

( ) i i f V x , is
estimated as the inverse of the lead time.

For a feasible solution, x i , we consider a weak neighbourhood and a heavy neighbourhood. The search is firstly carried out into the weak neighbourhood and then in the heavy one.

The weak neighbourhood consists of the interchanging of operations within the bottleneck station, E, and is named ( )

E i V x
. This weak search is carried out during

N w iterations.
The heavy neighbourhood consists of the interchanging of operations between any pair of stations and is named ( ) i V x . This search is carried out during N H iterations.

Every neighbour of any neighbourhood is depicted as ( )

i i V x .
In both cases the search for a better solution inside a neighbourhood is attempted during N loop iterations. 4 .

TABLE 4

The heuristic was initiated starting from the MH initial solution, from the solution provided by the MRPWH and from three other slight modifications from these two scheduling proposals.

We implemented 700 iterations for the weak neighbourhood search (N w ), 300 iterations for the heavy neighbourhood search (N H ) and 40 iterations were carried out inside the searching loop (N loop ), obtaining a computational time equal to approximately 950 seconds per each search process.

Table 5 shows the scheduling proposal according to the MSNSH heuristic method. Changes are mainly incorporated in stations 3 and 4 with respect to the proposal of the MRPWH method. These stations are critical and correspond to the bottleneck station identification during the computation of the neighbourhood search algorithm.

TABLE 5

Results show an improvement in the lead time equivalent to a reduction of 34 seconds (the new average lead time is assessed at 0:06:55).

Simulation and comparison of proposals using ARENA

We analysed the productivity and other performance indexes of the different proposals by developing a simulation model using ARENA simulation software, [START_REF] Kelton | Simulation with ARENA[END_REF]. Simulation software is especially useful for testing the efficiency of a new schedule prior to its implementation in real life, and allows the statistical data to be considered. Additionally we could contrast the level of accuracy of our simulation model by comparing it with the real situation that was implemented in the factory at that time.

We supposed the following hypothesis: the workers have all the material and required tools at their disposal, and the material is supplied to the line without interruption.

Times of operations were defined as having a normal distribution following Table 1. A whole labour turn of 8 hours was simulated. This time horizon allows for the variability of the results. The operations described in Table 1, together with the waiting queues that have to be considered due to possible delays (it is usual that a worker has to wait for the other worker in order to carry out a subsequent operation), are modelled using the constructive modules of ARENA.

The simulations of the MH initial scenario were particularly useful because this scenario allowed the validation of the simulation model (the results from the simulation were in accordance with the actual operation in the factory plant).

Simulation model

A simulation model was constructed for each specific solution. They all included some common parts and other different parts, depending on the specific scheduling of each solution. Figure 5 details the operation of the four stations of the assembly line for each case.

FIGURE 5

The first module defines the station, and the second module assigns the motorcycle in process to the station. The motorcycle then enters the corresponding module Station that is defined in the following figures 6 and 7, and which includes the operations that are carried out inside each station. Once the operations have been carried out, an "assign" module follows, stating that the station operations on the module where it waits for the remaining operations to be completed. The "route" module transports the entity from one station to the next.

The main building blocks are included in the Station module where the operations are implemented. Figure 6 and Figure 7 represent the simulation model for the schedule provided by the MSNSH heuristic.

FIGURE 6 FIGURE 7

Figure 7 addresses the operations in station 4, the most complex of the stations, and where the majority of operations are realized. The bottleneck of the production rate is described, as the following results show. The operations have been represented by "process" modules, where the resources are specified (which worker carried out the operations), and the operation time has a statistical normal variable.

Taking into account that workers can carry out operations at the same time, the entity is duplicated by a "separate" module named 'Duplicated'. Finally, both entities are rejoined by a "batch" module named 'Junction'. The figures include, together with the operation modules, some "hold" and "assign" modules that are included to represent the precedence relations and the waiting times between workers.

The processes for the MRPWH heuristic and the MH initial scenario are similar, apart from modifying the location of the operations that are interchanged in the scheduling, as can be seen in Tables 2, 3 and 4.

The simulation animation module is shown in Figure 8.

Comparison of results

This subsection shows the results for three schedules, focusing on the operation time for each station, the waiting times between workers and between stations, and the activity of the workers. The analysis of these three cases allows the efficiency and suitability of each schedule to be analysed.

Figure 9 depicts the operation time for each proposal analysed in each station.

MRPWH and MSNSH heuristics clearly outperform the initial proposal, providing a more balanced assembly line. In the initial situation the bottleneck was located in station 2, whereas MRPWH and MSNSH proposals locate the bottleneck in station 4, reducing the lead time drastically. However, some replications located the bottleneck in station 3 due to the variability (normal distribution) of the tasks' operation times.

By analysing Figure 9, we can also appreciate that MSNSH slightly improves the operation time provided by MRPWH.

FIGURE 9

It has to be taken into account that the location of the bottleneck is not particularly relevant, due to the individual nature of the assembly line conveyor, which moves in a discontinuous manner. The conveyor is consequently not activated until all the operations are completed, resulting in it being impossible to carry out work in advance, because the work in process (which is fastened to the line) cannot be displaced.

The idle times vary depending on the nature of the wait. Figure 10 shows the different waiting times of workers in the same station (top side), which occur when one worker has to wait for another worker to complete a task, as well as the waiting times that occur between stations (down side), appearing due to a non-exact balance.

FIGURE 11

Finally, the production level of motorcycles can be appreciated in Table 6 and Figure 12. Both MRPWH and MSNSH schedules produce a much larger number of motorcycles per day, with the MSNSH proposal producing the better values, with a higher average value, and a lower standard deviation, which is another important factor.

TABLE 6 FIGURE 12

Conclusion

This paper deals with a motorcycle factory real case study that corresponds to a discontinuous assembly line balancing problem. We implemented a traditional heuristic method that was adapted to the specific characteristics of a motorcycle assembly line (MRPWH). The heuristic proposal was compared with the output of a neighbourhood search method (MSNSH). We constructed a simulation model to test the different proposals. The simulations of the initial scenario permitted the validation of the simulation model.

The results from the simulation were in fact in accordance with the factory plant's current operation.

Both MRPWH and MSNSH provided much better solutions than the model implemented in the initial scenario. Also, the MSNSH proposal was slightly better than the MRPWH proposal, creating a greater productivity level and a more balanced line, which affected the operation times in the stations as well as the activity of the workers. Furthermore, the MSNSH schedule was so balanced than the bottleneck varied from station 3 to station 4, depending on the simulation replication, although it was mainly located in station 4. The measurement of the tasks' time by a statistical distribution allows a real representation of the current situation in the factory.

MH is at present implementing the MSNSH schedule in its factory, obtaining a productivity level of between 55 and 59 motorcycles (most of the days it is 56).

However, the efficiency of the workers is continuously improving, so the data in ( )

E i V x ( ) i i V x Else Random selection in ( ) i V x ( ) i i V x
While neighbour feasibility test fails Calculate ( )

( ) i i f V x If ( ) ( ) i i f V x < f(Best_V(xi)) Best_V(xi) ( ) i i V x While [f(Best_V(xi))> ( ) ( ) i i f V x ] OR Iterations < N loop If f(Best_V(xi))< ( ) ( ) i i f V x xi Best_V(

  been successfully implemented in the factory. Simulations are shown for all the cases, thus allowing the comparison between proposals.

  the whole vehicle; and the front and rear brakes to stop the motorcycle.The FURIA XP6 motorcycle assembly line is fed from manufacturing cells located near the corresponding station. Figure 2 depicts the layout of the manufacturing zone of the factory, where the assembly lines are located. Because Due to this model representing the largest part of the company's sales, there is one specific line dedicated to the model: assembly line 1. The other line is used for the company's remaining models. The assembly lines are operated in a discontinuous manner. The tasks are done with the chain stationary, and once the workers have finished their tasks the chain advances.

FIGURE 2 Figure 2

 22 FIGURE 2

  FIGURE 3

  Feasible solutions are encoded as Figure4depicts. The figure shows the encoding for the scheduling, corresponding to the factory's initial scenario. Each feasible schedule corresponds to a feasible solution, where three arrays are used for each possible worker distribution (right or left worker and both workers) at each station. The number of registers for the arrays of the stations must be sufficient, with the remaining non-used registers equal to blanks.

  Feasibility controls are required due to the precedence or position constraints. Only feasible solutions are considered, with the others being disregarded. The pseudocode of the algorithm is presented in Table

  completed in that station. Now the entity enters a "hold"

  FIGURE 8

  FIGURE 10

Figure 11

 11 Figure 11 represents the activity of the workers in the assembly line. The

FigureFigure 6 .

 6 Figure 1. FURIA XP6 motorcycle model

Figure 10 .

 10 Figure 10. Waiting times in the assembly line

TABLE 2

 2 As can be seen, the initial situation depicts an assembly line that is poorly balanced, with two stations (1 and 2) much busier than the others (3 and 4). The expected lead time is 0:12:49.

	Additionally,
	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

TABLE 3

 3

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

Table 1

 1

	F o r
	P
	e
	e r
	R
	e
	v i e
	w
	O n l
	y

is being revised. Additionally, new components are included in, or removed from, the motorcycle models each year. Due to all of these factors, the balancing of the assembly line must be continuously monitored and rebalanced in order to achieve better results and attain a state of continuous improvement.

Table 1 . Operations in the MH assembly line as they were carried out in the initial scenario Operation Precedence Position Average time Standard deviation

 1

	Operation	Precedence	Position Average time Standard

Table 2 . Initial scheduling in the assembly line

 2

	Worker

Table 3 . MRPWH scheduling proposal

 3

	Worker

Table 4 . MSNSH algorithm pseudocode

 4

	w=1
	Do
	Do
	Do
	Do
	If w=1	Random selection in

Table 5 . MSNSH scheduling proposal

 5

	xi)

Table 6 . Daily productivity results

 6

		Daily average productivity	Standard deviation	Interval
	MH initial scenario	35.20	0.51	(34,36)
	MRPWH heuristic	56.25	0.70	(55,57)
	MSNSH heuristic	57.05	0.44	(56,58)

http://mc.manuscriptcentral.com/tprs Email: ijpr@lboro.ac.ukInternational Journal of Production Research