Jesus ; Gonzalez-Feliu
email: jesus.gonzales-feliu@let.ish-lyon.cnrs.fr

Jean-Louis Routhier
email: jean-louis.routhier@let.ish-lyon.cnrs.fr

Charles Raux
email: charles.raux@let.ish-lyon.cnrs.fr

An attractiveness-based model for shopping trips in urban areas

Keywords: urban goods movement, shopping trips, simulation, catchment area model, attractiveness

This paper presents a modelling approach to characterise private car shopping trips in a city logistics point of view, in order to connect these movements with those belonging to urban freight distribution in the supply chain. The proposed modelling framework is a two-step procedure articulated as follows. First, an attraction model estimates the number of private car shopping trips arriving to each section of a given urban area. Second, a catchment area model relates the shopping trip destinations with the household locations. The model is calibrated using the data of both a database of the commercial activities and the recent household trip survey made in 2006 in Lyon urban area (France). We present the main results issued of the various simulations as well as several application examples and proposals, most of them made in a public policy perspective.

INTRODUCTION

Shopping trips are an important component of urban traffic, since they are related to movements of people and goods: in France they represent about 50% of the total PCU km travelled by urban goods (including commercial vehicles). Moreover, travel patterns for shopping are subjected to important mutations, mainly because of the urban developments and trends on new retail and logistic distribution schemes.

Shopping trip estimation is usually made related to people, using the classical 4-steps modelling approach: it follows the hypothesis that these trips depend on the same factors as other trip purposes. However, several studies show that shopping trip generation rates are related to factors that are different from those of work trips for instance. Moreover, shopping trips are difficult to characterise due to the complexity and the variety of trip chains where at least one shopping activity takes place, and the available data sources have not been created for specific shopping trip generation.

We propose a new modelling approach that relates shopping trip flows to retailing supply. To do this, we calculate the catchment area of each zone i, in terms of number of households of each zone j that will make a journey that contains at least one shopping stop at i. The approach uses standard household travel survey data and available data about retail supply.

First, we present the main background issues to settle the context of our research. Then, the model is presented, focusing respectively on the definitions, assumptions and other hypothesis taken into account, the attraction phase and the catchment area phase. After presenting the two phases of the model, we present the calibration procedures and their main results. Finally, two examples of applications are proposed. More precisely, we show a simulation of extreme location and distribution scenarios in the urban area of Lyon (France) and an analysis of the impacts of e-commerce and teleshopping practices, in progressive usage increasing scenarios. The scenarios are compared on the criteria of total PCU km travelled by urban goods and total emissions of CO 2 , taking into account the different types of vehicles used, i.e. private cars, light city freighters, urban trucks and semi-articulated vehicles.

BACKGROUND ISSUES

In the last decades, city logistics has been developed to deal with the main problems of urban freight distribution, studying freight movements in urban areas and proposing solutions to reduce congestion and pollution as main problematics. Recent studies have defined and characterised the different movements of urban goods (Patier, 2002;[START_REF] Ségalou | The environmental assessment of urban goods movement[END_REF][START_REF] Russo | Demand models for city logistics : a state of the art and a proposed integrated system[END_REF]. Two categories are predominant and represent about 90% of the overall urban goods movements: establishment supply movements, which are related to freight distribution between the different activities, and end-consumer commodity movements, where the purchased goods are moved by the consumer, related to shopping trips. The first group of urban freight movements, which corresponds to the exchange of goods between different establishments, has been one of the most studied subjects in city logistics research.

End-consumer movements, less studied in urban freight transport science, is an important component of urban freight transport, representing about 50% of the total km.PCUs (Passenger Car Units) of urban freight transportation, where establishment supply movements represent only 40 % (Patier, 2002). The remaining 10% contains the city maintenance and construction logistics movements, the waste distribution and a number of small particular activities such as postal services, among others In France, shopping trips represent 10% of the total person trips in working days, i.e. from Monday to Friday, and 25% in Saturday [START_REF] Routhier | Mesurer l'impact du transport de marchandises en ville -le modèle Freturb (version 1)[END_REF]. In terms of pollution, the urban transport of goods, including end-consumer movements, produces about 25% of the total CO2 emissions for transport, 35% of the NOx emissions and 40 to 50% of the particulate [START_REF] Ségalou | The environmental assessment of urban goods movement[END_REF]).

3 Moreover, mutations related to the development of new technologies and consumption trends have an important impact on the shopping behaviour, which influences the different travel habits, so the characteristics of the shopping trips [START_REF] Routhier | ETHEL II -Rapport final[END_REF]. Land use policies developed by public administrators will have also an influence on commercial activities location and the services they propose [START_REF] Routhier | Mesurer l'impact du transport de marchandises en ville -le modèle Freturb (version 1)[END_REF]. Moreover, the location of retail activity areas has an influence on freight flow generation, in both sides of the extended supply chain (traditional freight distribution and end-consumer movements). Although the contribution of end-consumer movements to the urban flows of motorised vehicles is important, they are rarely taken into account in city logistics, at both policy making support and transport planning and optimisation levels.

Shopping trip estimation is generally made by classical four-step models [START_REF] Ortuzar | Modelling Transport[END_REF]Hensher and Button, 2001) but in other fields non transportation-based models can also be found (Till and Timmermans, 1992; Bawa abd Ghioshi, 1999; Long-Lee and Pace, 2005; [START_REF] Kubis | Analysis of Location of Large-area Shopping Centres. A Probabilistic Gravity Model for the Halle-Leipzig Area[END_REF].

Trip generation models relate trip generation rates to land-use and household characteristics. In general, the focus of research is the number of trips generated and their geographic distribution. Most of these models are integrated in other methods like the well-known family of classical four-steps models [START_REF] Ortuzar | Modelling Transport[END_REF], which is commonly used in general trip characterisation. The socio-economic characteristics of the trip makers are assumed to be significant determinants of travel behaviour because the factors determining the number of all types of trips are assumed to be the same as those for shopping trips [START_REF] Cubukcu | Factors Affecting Shopping Trip Generation Rates in Metropolitan Areas[END_REF]. These factors are, among others, income range, age, gender, employment status, car ownership, and household size. The physical and demographic characteristics of the area include: employment, population, and density. However, it is reasonable to believe that there are metropolitan area characteristics and trip maker characteristics which impact is significant only for shopping trip generation rates [START_REF] Cubukcu | Factors Affecting Shopping Trip Generation Rates in Metropolitan Areas[END_REF]. However, most of the studies dealing with shopping trips specifically are in general related to the regional or national level (Vickerman and Bramby, 1985;[START_REF] Badoe | Urban and Travel Changes in the Greater Toronto Area and the Transferability of Trip-generation Models[END_REF]Cubucku, 2001;Simma et al., 2005); few models are proposed to characterise shopping trip generation factors in a urban context (Ségalou, 1999;Black et al., 2007). Most of the works deal with econometric and empirical approaches based on surveyed data, mainly collected to calibrate this type of models.

Shopping trip distribution is more heterogeneous. The main categories are the classical trip distribution methods [START_REF] Ortuzar | Modelling Transport[END_REF] based on gravitary and entropyminimization approaches, and discret choice approaches (Thill and Timmermans, 1992). Moreover, in real estate research, catchment area models can also be used to estimate the main customer's locations in order to estimate trip distances and household (or inhabitant) travel behaviour categories (Long-Lee and Pace, 2005; Kubis and Hartman, 2007). However, these models are built on specific surveys for one store or one retailing area and have not been conceived to be used with standard data on several urban contexts. • An inbound trip, arriving at the considered commercial area. This is the part of the trip chain that is defined as shopping trip in classical passenger trip characterisation.

THE PROPOSED MODEL Definitions and model description

• An outbound trip, starting at the considered commercial area, where the goods that have been purchased are also travelling. Although the end-consumer goods movements take place only when the outbound trip starts, they are closely related to inbound trips, as both types of shopping trips are always made consecutively one after the other (see Figure 1). A particular type of shopping-related trips is that of trips which purpose activity at the origin and the destination is related to shopping. The Shopping-Shopping trips can be considered as inbound trips for the arriving zone and outbound trips for the departing zone.

Trip generation (attraction)

In this phase the shopping trips attracted by each zone are generated. To do this, we built a model using the general trip rate generation modelling framework.The model will be able to generate the overall shopping trip rates (all modes) at each zone, for both emission and attraction. In a first time, we will characterise the shopping trips at the shopping destination, which represents both the inbound trip's destination and the outbound trip's origin. We will note the number of private car trips at a shopping destination i the private car shopping attractiveness of i and noted i T .

Trips are the basic unit for zone-aggregation trip generation models, and the output variable is in general the number of trips which departure or arrival is located in a pre-identified zone i during a given time unit (Hobbs, 1979). As already said, the proposed model estimates the private car shopping trip attraction rates [START_REF] Ortuzar | Modelling Transport[END_REF]. Following the general trip generation rate modelling framework (Hobbs, 1979;Cubucku, 2001;[START_REF] Ortuzar | Modelling Transport[END_REF], the number of shopping trips T made by private car starting at or arriving to zone i can be formulated in the following way:

T i = f (Ret i , X i , Tech i )
where T i is the total number of shopping trips generated in section i; Ret i the set of retailing activity characteristics vector in section i; X i the set of socio-economic characteristics of people and households belonging to section i; and Tech i the set of technology characteristics. This equation can be formulated as a vector function: Other interesting databases are obtained from the household trip surveys, both local or national, or the Commerce and Industry Chamber national or regional surveys, although these data are less accessible and not available for all the medium urban areas.

T = f (Ret, X, Tech) (1) with                     = Z i T T T T T 2 1 ;                     = Z i Ret Ret Ret Ret Ret 2 1 ;                     = Z i X X X X X 2 1 ;                     = Z i Tech
The technology vector Tech is included in the framework only in a conceptual point of view. Current surveys do not have enough elements to well define tele-shopping usage trends, and other behavioural trends related to home delivery. Moreover, current data sources do not allow to characterise the relations between teleshopping trends and car usage behaviour (Beauvais, 2005). For these reasons we will not include in a first time the effects of teleshopping and home delivery services in our model, but only include them in a conceptual modelling framework.

The trip attraction function can be defined in different ways. We propose a multilinear function, defined as a pondered sum of the considered variables, and can be obtained by linear regression techniques applied directly on the available data. We can then define this function as follows:

∑ ∑ ∑ + + = e g g g f f e e e
Tech a X a Ret a T

We will estimate the coefficients by linear regression techniques, more precisely by the square minimum method.<Q Classical approaches are applied to estimate the total number of shopping trips in a urban area (Vickerman and Bramby, 1984;[START_REF] Cubukcu | Factors Affecting Shopping Trip Generation Rates in Metropolitan Areas[END_REF] or those trips for each zone of a urban area (Ségalou, 1999; UK grocery trips). A different approach (Gonzalez-Feliu et al., 2010) proposes a modelling approach on three categories of urban space:

• The main central urban area contains the main city of the urban region and sometimes also other urban suburbs which can be assimilated to the main city, because of a continuity of the urban landscape.

• The near periphery includes the urban zones limiting with the central urban area first ring.

• The rest of towns of the urban area belong to the far periphery.

The main problem of using this approach is that splitting the statistical population in three categories decreases the number of zones used to calibrate each of the three models, making these three sets of data not always enough accurate to ensure a quality econometrics analysis in the calibration phase of these models. For these reasons, we will test both approaches and compare them.

Catchment area distribution

Once the inbound shopping trips are generated at their destinations, we need to define the origin of these trips and the destination of the related outbound trip. From the household trip surveys results, we observe that households correspond to 80% of the inbound origins and to 85% of the outbound destinations. Therefore 70% of the shopping trips are included in household-shopping-household chains.

We propose to connect shopping destinations to households using a catchment area model. The proposed approach is a gravity probabilistic model. This model has a similar form than classic gravity distribution models (Ortuzar and Willumsen, 2001):

β hi i h i hi c T NrH A N = where ∑ = k hk h i c NrH A β 1
In order to take into account the retailing supply of both household-related and destination zones, we introduce two other variables, one for each zone. The model can then be formulated as follows:

β α α hi i h h i i hi c T NrH E E A N 2 1 = where ∑ = k hk h h i i c NrH E E A β α α 2 1 1

MODEL CALIBRATION

The proposed framework is built and calibrated with the available standard data of the urban area of Lyon, which has about 2.000.000 inhabitants and 800.000 households. The main data sources are an extract of the register file of companies (SIRENE file) of the chosen area, the corresponding census database (INSEE file), and the 2006 personal trip survey, which follows a French standard [START_REF] Certu | L'enquête ménages déplacements 'standard CERTU' : Guide méthodologique[END_REF]. In this survey, the urban community of Lyon is divided into several small zones, grouped into macro zones. However, it is important that the methodology will be able to be applied to a chosen urban area or generalised into a model able to be used on any urban area without the need of a calibration and a testing phase. Shopping trip behaviours are heterogeneous and difficult to characterise using microsimulation approaches, as the available data from surveys has not been collected for shopping trip characterisation and the high cost of these surveys made them adaptable to develop macro-simulation models, but not to micro-simulation ones. The model parameters are calibrated separately for each phase, following a mixed approach that combines an econometric analysis (using linear regression techniques) and the fitting method proposed by [START_REF] Hyman | The calibration of trip distribution models[END_REF]. In the following subsections we describe the calibration methods and the main results of the calibration analysis.

Attraction phase

The attraction phase uses a multilinear function. To calibrate it, we define each coefficient by linear regression. We have implemented the model using R Commander [START_REF] Fox | The R Commander: A Basic-Statistics Graphical User Interface to R[END_REF], combining the different variables in order to obtain the most precise model. The model can be formulated as follows:

1C Model:

j j j j
+ - + - + + + =
We observe that the R² is close to 1 in all three categories of urban space. Moreover, the F test is close to 0 in two of the cases (see Erreur! Source du renvoi introuvable.). However, the one-category model is more robust, due to the highest number of individuals considered in the linear regression procedure.

Catchment area phase

The catchment area phase presents various parameters that must be calibrated. We propose a calibration method that combines an econometric analysis (Kubis and Hartman, 2007) with a fitting procedure based on that of [START_REF] Hyman | The calibration of trip distribution models[END_REF]. This calibration method is presented as a sequential procedure that can be resumed as follows:

1. Determination of the coefficients for both small and big stores by linear regression.

2. Revision of the model and final formulation 3. Determination of the coefficients to calibrate using the fitting procedure and iniziatilation of the procedure 4. Iterations for coefficient improvement using the fitting method

The linear regression analysis starts from the following formulation for the catchment area model:

β γ α α hi i h h i i hi c T NrH E E A N 2 1 =
We approximate in a first time by excluding in the analysis the A ij elements. Then, we can linearize the function using the logarithmic operator:

hi h h i i c NrH E E K T P log log log log ) ( log 2 1 β γ α α + + - + =
Finally, using R Commander [START_REF] Fox | The R Commander: A Basic-Statistics Graphical User Interface to R[END_REF] we establish the form of the model as well as the variables and coefficients that are included in the final formulation. We found that the two best approximations are the following: The commercial supply for small stores in the residence zone can be supposed as a negligible variable, as in the different combinations that have been tested the coefficients were close to 0 and less influent than the commercial supply for small stores in the shopping zone. The number of households has an inverse effect as that we can intuitively hypothesise, but this can be seen in the following way: the influence of household is pondered by the distance. In fact, people trends are to go shopping close to the household zone when the commercial supply satisfies their needs, and the intra-zone distance is in general smaller than the inter-zone distance.

Then, the final model can be formulated. Taking into account the similarity of the distance coefficient, we propose the following catchment area model:

β γ α α γ α α hi b h big b h b i big hi s h small s h s i small hi hi c NrH K E E A NrH K E E A T         + = 2 1 2 1 with ∑ ≠ = i k hk s h small s h s i small i c NrH K E E A β γ α α 2 1 and ∑ ≠ = i k hk b h small b h b k big hi c NrH K E E A β γ α α 2 1
The third step is the fitting procedure. We start the procedure by initialising the algorithm. To do this, we first select the parameters that will be recalibrated. Alpha and gamma belong to the A ij elements, and they should be calibrated using other approaches, such as that proposed by [START_REF] Ortuzar | Modelling Transport[END_REF], but they need complete O/D matrices where the total number of emission trips and the total number of attraction trips are the same quantity, which is not our case. For this reason, we keep the coefficients estimated by linear regression and make the calibration only on beta.

We define the following mean values of c [START_REF] Hyman | The calibration of trip distribution models[END_REF]: a. We make m=m+1. Using the latest value for β, noted β m-1 , we calculate the catchment area matrix using the model ( 3) and obtain the new mean modelled trip distance cost cm-1 using the following expression:

∑ ∑ = hi R hi hi ij R hi
∑ ∑ = hi hi hi hi hi N c N c ) ( ) ( ) ( β β β We calculate ) , ( 1 R ij m ij N N - δ
. If this value is sufficiently close to accept the precision of the model, the iterations are stopped, otherwise we go to step b. b. We improve β in the following way:

1 1 1 1 ) * ( ) * ( - - - + - - - - = m m m m m m m c c c c c c β β β
The iteration phase is repeated until the precision of the model is accepted as good, i.e. when cm is supposed to be sufficiently close to c*. Moreover, if the parameter is not enough improved after I iterations, we also stop the calibration process. We report in Table 3 the main calibration results of the model. After 50 iterations, we obtain an error on the mean distance lower than 5%. We observe that in a global perspective, the results are quite satisfactory. The distribution results are less performing but the estimation of the total distance due to private cars shopping trips can be estimated with a small error.

APPLICATION DOMAIN AN USAGE EXAMPLES

Possible uses of this can be related to both passenger and freight transportation, and also to marketing and real estate fields. For a public policy usage, different choices and applications can be proposed. For example, the proposed model can be used to simulate strategic planning scenarios to estimate the impact of urbanistics, legislation and other policy actions in both people and the global freight movement trends.

Another use can be to study the limits of several policies (Gonzalez-Feliu et al., 2010). It is important to note that the new distribution services (home delivery, proximity reception points, etc.) have a direct impact on shopping trips. In order to simulate the substitution of current practices by these new distribution services, the shopping trips by private car have to be characterised in order to define and simulate the substitution trends (from current shopping trips to the new distribution routes). [START_REF] Routhier | ETHEL II -Rapport final[END_REF] propose four extreme urban retailing supply scenarios:

1. concentrated shopping malls; 2. traditional medium level stores; 3. e-commerce plus home delivery; 4. e-commerce plus local pick up points.

To simulate, the inter-establishment movements (IEM) have been estimated using the LUTI model Freturb (Routhier and Toilier, 2007) and the end-consumer movements (ECM) using the proposed shopping trip model. The substitution of shopping trips for the teleshopping scenarios is made using the procedure proposed by [START_REF] Routhier | ETHEL II -Rapport final[END_REF]. Moreover, the CO2 emissions can also be estimated from the travelled distances and an estimation of the mean speed in each section using the commercial software Impact Ademe (ADEME, 2003). The "all hypermarkets" situation is a little better than the reference, in terms of freight flows for both types of movements. The establishment supply ones are drastically reduced because the logistics systems of the commercial centres allow a better rationalisation of freight distribution flows (about -87% of the km.PCU respect to the reference situation). However, the end-consumer movements do not increase with respect to the reference(we observe a light decrease, near -3%), due to the need of using the private car to arrive to the peripheral commercial centres that a priori increases the length of private car shopping trips and the concentration of the activities in this form of commercial supply that reduces the number of trips. The global result is a reduction of a little more than 15% of the total km.PCU, which shows the importance of end-consumer movements in the overall logistics system of a urban community. The "all proximity stores" scenario has a higher impact. Although the establishment supply movements increase (+90% of the km.PCU), the end-consumer movements decrease because of the proximity that incentives the "on foot" shopping trips (about -87%). This is traduced by a reduction of the total km.PCU of 56% respect to the reference. The "all home delivery services" scenario has an impact on the end-consumer movements: although no shopping trips are observed, considering everybody is delivered home, the home-delivery trips (another component of end-consumer movements) present an important increasing trend, which is traduced by a reduction of -85% in the end-consumer km.PCU. The distribution to the logistic platforms for home delivery has been considered to be made in a similar way respect to the hypermarket distribution, in a very rationalised way (a collaborative situation is simulated), so the establishment supply movements are close to those obtained for the first scenario. The same hypothesis are made to the fourth scenario, but the reception points system has a better optimised distribution to the reception points. Moreover, we obtain a small number of private car shopping trips. This is derived from the current proximity shopping behaviour, considering than several households will use the car to take big quantities of goods from the reception points. The end-consumer movements in this situation have a similar trend respect to the "all small stores" one (-93%). These trips contain both the shopping and the reception points delivery movements. The overall freight flow trends for these two scenarios are respectively -86% and -92% respect to the reference situation. We see than the home delivery services are not the best solution, since the transportation system that replaces shopping trips is very rigid and its efficiency improvement is difficult (Alligier, 2007). The "all reception points" scenario seems to have a more positive impact, because a big reduction of the travelled distances can be made for both types of movements, which is not the case in the first family of scenarios. A combination of these four can be used by stakeholders to define urban commercial supply policies that integrate both types of movements in a city logistics point of view.

The catchment area approach can also be useful to simulate e-commerce policies and behavioural issues related to the new distribution forms. Although specific data about ecommerce behaviour in the zone of simulation is difficult to found, the findings of Rohm and Swaminathan (2004) can be applied to medium and big urban areas. These authors define 4 types of online shoppers on the basis of a specific survey for e-commerce behaviour on standard categories of population. We can state than although the store-oriented shoppers are only 15% of the sample, a catchment area approach including variables such as the commercial supply and the geographical distance (among others) represent the vehaviour of 3 categories, i.e. 89% of the population. For this reason, assuming the usage rates of ecommerce, the proposed model can distribute the related trips in order to reconstitute the delivery routes for several e-logistics schemas [START_REF] Routhier | ETHEL II -Rapport final[END_REF][START_REF] Durand | Pas de developpement durable du B to C sans veritable e-logistique urbaine[END_REF].

We propose several results for studying the effects of e-commerce introduction. The usagerate of these services is increased from 0% (assuming than in 2006 the e-commerce usage rates are negligible) to 100% with an increase of 10% for each scenario. The attractiveness model is used to estimate the potential demand of e-commerce distribution services for each zone and the remaining classical shopping trips. Following the simulation methodology of [START_REF] Durand | Pas de developpement durable du B to C sans veritable e-logistique urbaine[END_REF] for route simulation, we can estimate the travelled distances of these forms of distribution. Then, the total distances (e-commerce distribution and classical shopping trips). The simulated data are the following: We observe that e-commerce competition scenarios are less performant that their analogous collaborative scenarios. Indeed, with an "all home delivery" scenario the total reduction of the number of km. PCU is near -60% whereas in a collaborative situation it is about 85%. The "all reception points" scenario presents similar results (-65% with respect to -92% for the collaborative scenario). We also observe than in a situation without important behavioural changes, a 50% reduction of the total number of Km.PCU is reached only with a 80% of usage rates for e-commerce.

CONCLUSION

In this paper we presented a new modelling approach for shopping trip simulation, in order to reproduce an estimation of end-consumer goods movements which would be able to be compared to freight distribution trips. Considering only the private car shopping trips, which are those that interact with freight distribution trips, we propose a two-step model. First, the attracted shopping trips are generated at each shopping destination section. Second, each shopping destination section is related to the possible household locations of the shoppers using a gravity catchment area model.

Another important element of shopping trips in urban areas is that each zone does not have the same characteristics. Following the administrative typology of rings to classify the urban community's towns, we divided the considered zones into three categories for trip generation: the Central Urban Area (CUA), which is the main city's area, the Near Periphery (NP), which represents the first ring, and the Far Periphery (FP), which contains all the other towns of the urban community. Similarly, we also propose a two-category model for the catchment area phase, in order to take into account the differences between shopping trips made for purchasing at small retailers and those related to big stores and hypermarkets.

The shopping trips at each retailing zone is well estimated, and similar results are obtained for the outbound trips. However, the inbound trips are less well estimated due to the influence of working activities (no information about the overall number of employees for all the economical activities is included). Moreover, two examples of modelling applications made with the proposed model are presented.

In order to improve the model, and to define an approach adaptable to each city, we have applied this methodology to a data file containing zones of two urban areas of different size.

The main results are encouraging (similar regression coefficients and good approximations in both urban areas), and a deeper study is in process. Moreover, a better characterisation of the Commercial Centres, as for example defining them by their surface or their total revenue, have to be studied.

  The proposed model follows the definitions and principles enounced in Gonzalez-Feliu et al. (2010). In this work; the authors propose to model directly the private car shopping trips in order to compare them to establishment supply flows. According to several researches (Ségalou, 1999a; Dablanc et Pecheur, 2000; Toilier et al., 2005; Chambre de Commerce et d'Industrie de Lyon, 2008; Gonzalez-Feliu et al., 2010), we can observe than shopping trips are included in more complex trip chains. In all these chains, commercial activities are related to two connected shopping trips:

Figure 1 :

 1 Figure 1: schema of the shopping-related trips

2 1

 2 Each vector contains a set of possible variables that can be included in the shopping trip generation model.For this research the variables considered in the vector Ret i are the following: NrSMC : Number of small and medium stores; NrBS : Number of supermarkets and big commercial surfaces; NrVBS : Number of hypermarkets and similar stores; Emp-SMC : Number of employees in the section's small and medium stores; Emp-BS : Number of employees in the section's supermarkets and big commercial surfaces; Emp-VBS : Number of employees in the section's hypermarkets and similar stores; CC : Presence of an commercial centre(Ségalou, 1999), presented as a binary variable which takes the value 1 if at least one extra-urban commercial centre is located inside the section's perimeter and 0 otherwise. These variables are obtained by aggregation of the SIRENE file data, for the considered city. The SIRENE databases is a set of French establishments information data files, which 6 contains among others the retail activities' basic information (location, size category, commerce category, number of employees, etc.) for each city. The population socio-economic characteristics vector X includes: POP: Population of the considered section NrH: Number of households of the considered section DPOP: Population density; DH: Household density; These variables are extracted from INSEE 1 population census files. In France, each local administration has access to the population and household data of the national census.

Figure 2 :

 2 Figure 2: Maps of the Lyon urban area and the considered zones (source: L. Bouzouina, LET)

10

 10

N

  is the number of measured trips arriving at i which households are located in h.The fitting procedure can the be defined, and is composed by several stages that are defined as follows: 0. Initialisation: Starting on 0 β , a a better value of β is estimated in the following way:

Table 1 -

 1 Linear regression results for the best attraction analysis

	Model	R²	F Test	Individuals
	1C	0,82 1,28.10-9	33
	3C -CUA	0,86	0,31	8
	3C -NP	0,77	0,007	14
	3C -FP	0,95	0,002	11

Table 2 -

 2 Linear regression results for the best catchment area model

	Model	R²	F Test	Individuals
	Small retailing activities		0,47 1,04.10 -23	179
	Big stores		0,49 3,98.10 -25	176

Table 3 -

 3 Estimation results on the data ofLyon (2006)

	Surveyed	Estimated Error

Table 4 -

 4 Simulation results for progressive increase of e-commerce user rates

	Scénario		Aval Livraison Aval Achats Total
	0		0	25955939	25955939
	LAD	10%	4%	-10%	-6%
		20%	8%	-20%	-12%
		30%	13%	-31%	-18%
		40%	17%	-42%	-25%
		50%	21%	-52%	-31%
		60%	25%	-63%	-38%
		70%	29%	-74%	-44%
		80%	33%	-84%	-51%
		90%	38%	-95%	-57%
		100%	40%	-100%	-60%
	PR	10%	2%	-8%	-6%
		20%	4%	-17%	-13%
		30%	6%	-26%	-20%
		40%	7%	-34%	-27%
		50%	9%	-43%	-34%
		60%	11%	-52%	-41%
		70%	13%	-61%	-48%
		80%	15%	-70%	-55%
		90%	17%	-79%	-62%
		100%	18%	-82%	-65%

th WCTR,July 11-15, 2010 -Lisbon, Portugal

ACKNOWLEDGEMENTS

Part of the proposed methodology was developed in the context of ETHEL II project, part the GICC, a research program on climatic change financed by the French Agency of the Environment and the Energy (ADEME).

The authors should also like to acknowledge Florence Toilier, studies engineer at ENTPE, Vaulx-en-Velin (France) for her help in the data collection and the estimation of freight distribution trips using FRETURB, as well as Frédéric Henriot, PhD. student at Université de Lyon 2 for his help and exchanges in the constitution of scenarios.


Ph.D Liesbeth Ceelen
email: lceelen@vub.ac.be

Ward De Spiegelaere

Michael David

Jurgen De Craene

Mathieu Vinken

Tamara Vanhaecke

Vera Rogiers

Critical selection of reliable reference genes for gene expression study in the HepaRG cell line

Keywords: HepaRG, RT-qPCR, reference gene, normalization

The human HepaRG cell line has shown to be a valuable in vitro tool for repeated exposure to chemical compounds and to evaluate their potential toxic outcome. Seen the importance given by the actual EU legislation of cosmetics and chemical substances to the use of in vitro methods in human safety evaluation, one can expect that HepaRG cells will gain importance as human-relevant cell source. At the transcriptional level, RT-qPCR assays are often used to obtain quantitative results. The choice of internal control is important since it may affect the study outcome. Indeed, it is well-known that expression levels of traditional reference genes can vary across tissue types and across experimental settings within one specific tissue type.

From a review of the scientific literature, it appears that, for HepaRG cells, S18 often is used as internal control, but without any evidence of its expression stability in this cell line.

Therefore, we aimed to select the most optimal reference genes for gene expression studies in HepaRG cells and to check whether S18 is a suitable reference gene. Twelve candidate genes' expression stability level was analyzed by three algorithms (geNorm, BestKeeper, Normfinder), which identified the optimal single reference gene (TBP) and the most suitable set of reference genes (TBP, UBC, SDHA, RLP13, YHWAZ, HMBS, B2M and HPRT1) for HepaRG transcriptional profiling. This study provides a new set of reference genes that is suitable for testing whenever RT-qPCR data for HepaRG cells are generated. The most stable ones can then be selected for further normalization.

A c c e p t e d M a n u s c r i p t

Introduction

Drug-induced liver toxicity is among the main reasons to withdraw drugs from the market [START_REF] Lasser | Timing of new black box warnings and withdrawals for prescription medications[END_REF].

This may be partially explained by the poor capability of the currently available animal-based in vitro systems to predict human toxicity in preclinical testing and hence implicates an urgent need for a reliable in vitro system that appropriately reflects the human in vivo situation for hepatotoxicity testing. Notwithstanding primary human hepatocytes still represent the gold standard in the field of liver-based in vitro models, many efforts have been done for extended periods of time, to develop a human hepatic cell line expressing various liver-specific functions at an in vivo-like level. The hepaRG cell line, derived from a human hepatocellular carcinoma, seems to meet the majority of these criteria and therefore represents for the time being a suitable alternative for freshly isolated human hepatocytes in drug screening [START_REF] Gripon | Infection of a human hepatoma cell line by hepatitis B virus[END_REF][START_REF] Guillouzo | The human hepatoma HepaRG cells: a highly differentiated model for studies of liver metabolism and toxicity of xenobiotics[END_REF][START_REF] Kanebratt | Evaluation of HepaRG cells as an in vitro model for human drug metabolism studies[END_REF].

Seen the importance given by the actual EU legislation of cosmetics (76/768/EEC and 2003/15/EC) and chemical substances to the use of in vitro methods in human safety evaluation, it could be expected that in the near future the HepaRG cells will gain importance as human-relevant cell source. Actually, the human HepaRG CYP induction test method is currently under validation at the European Centre for the Validation of Alternative Methods (ECVAM Technical report on the Status of Alternative Methods for Cosmetics Testing, 2008Testing, -2009; http://ecvam.jrc.ec.europa.eu/). The HepaRG cell line has indeed shown high metabolic capacity compared to primary human hepatocytes [START_REF] Guguen-Guillouzo | Stem cell-derived hepatocytes and their use in toxicology[END_REF][START_REF] Guillouzo | Evolving concepts in liver tissue modeling and implications for in vitro toxicology[END_REF]. These cells have been used for longterm repeated exposure to evaluate the potential toxic effects of different chemical entities [START_REF] Guguen-Guillouzo | Stem cell-derived hepatocytes and their use in toxicology[END_REF][START_REF] Josse | Long-term functional stability of human HepaRG hepatocytes and use for chronic toxicity and genotoxicity studies[END_REF]. Once differentiated using 2% dimethylsulfoxide (DMSO), cells at confluency show indeed a relative preserved expression of transcripts encoding various phase I, phase II and antioxidant enzymes, membrane (apical, canalicular and basolateral) transporters, the nuclear receptors constitutive androstane receptor and pregnane X receptor, aldolase B and albumin, up to 6 weeks in culture [3-5, 7, 8]. Expression of a plethora of markers/functions has been A c c e p t e d M a n u s c r i p t evaluated both at the transcriptional and the translational level. With respect to the former, reverse transcription quantitative PCR (RT-qPCR) assays were generally adopted often using the S18 rRNA gene as single reference gene [START_REF] Josse | Long-term functional stability of human HepaRG hepatocytes and use for chronic toxicity and genotoxicity studies[END_REF][START_REF] Vee | Functional expression of sinusoidal and canalicular hepatic drug transporters in the differentiated human hepatoma HepaRG cell line[END_REF][START_REF] Aninat | Expression of cytochromes P450, conjugating enzymes and nuclear receptors in human hepatoma HepaRG cells[END_REF][START_REF] Antherieu | Stable expression, activity, and inducibility of cytochromes P450 in differentiated HepaRG cells[END_REF][START_REF] Cerec | Transdifferentiation of hepatocyte-like cells from the human hepatoma HepaRG cell line through bipotent progenitor[END_REF][START_REF] Dumont | Preferential induction of the AhR gene battery in HepaRG cells after a single or repeated exposure to heterocyclic aromatic amines[END_REF][START_REF] Lambert | Reproducible chemical-induced changes in gene expression profiles in human hepatoma HepaRG cells under various experimental conditions[END_REF].

In the entire qPCR process, a number of variables, including RNA extraction (yield, integrity and DNA contamination), efficiency of reverse transcription and PCR cycling impose the necessity to include a number of controls to guarantee the integrity of every step and hence, to achieve reliable and accurate results. It is therefore important to use reliable reference genes as internal controls in order to normalize mRNA data [START_REF] Teste | Validation of reference genes for quantitative expression analysis by real-time RT-PCR in Saccharomyces cerevisiae[END_REF][START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF][START_REF] Thellin | Housekeeping genes as internal standards: use and limits[END_REF][START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF]. Their efficacy, however, must be experimentally validated for a particular tissue type and/or experimental condition. Indeed, it is known that the expression profiles of commonly used control genes can vary across different cell types, but also within one cell type subjected to a different experimental treatment [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF][START_REF] Thellin | Housekeeping genes as internal standards: use and limits[END_REF][START_REF] Nolan | Quantification of mRNA using real-time RT-PCR[END_REF][START_REF] Gonzalez-Verdejo | Selection of housekeeping genes for normalization by real-time RT-PCR: analysis of Or-MYB1 gene expression in Orobanche ramosa development[END_REF][START_REF] Spiegelaere | Quantitative mRNA expression analysis in kidney glomeruli using microdissection techniques[END_REF]. Consequently, the selection of inappropriate reference genes may significantly affect the study outcome by inaccurate data interpretation as highlighted in a number of previously reported studies involving a variety of cell types and experimental situations [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Ferguson | Impact of Reference Gene Selection for Target Gene Normalization on Experimental Outcome Using Real-Time qRT-PCR in Adipocytes[END_REF][START_REF] Nicot | Housekeeping gene selection for real-time RT-PCR normalization in potato during biotic and abiotic stress[END_REF].

Therefore, the current study was set up to search for the best reference genes for RT-qPCR assays in the HepaRG cell line using the geNorm [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF], BestKeeper [START_REF] Pfaffl | Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations[END_REF] and Normfinder [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF] software. The quality assessments described, are in line with the recently published MIQE guidelines (Minimum Information for Publications of Quantitative Real-Time PCR Experiments) that provide authors and reviewers with the minimum required information to ensure the quality of interpretation and repeatability of qPCR experiments [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF].

Materials and methods

2.1.HepaRG Cell cultures

HepaRG cells, originally derived from a liver tumor of a female patient suffering from hepatocarcinoma [START_REF] Gripon | Infection of a human hepatoma cell line by hepatitis B virus[END_REF], were purchased from Biopredic International (BPI, Rennes, France). For A c c e p t e d M a n u s c r i p t the present study, HepaRG cells were cultured as previously described. Thus, cells at three different passages were cultured at low density (i.e. 2.6 × 10 4 cells/cm 2 ) or high density (i.e. 4.5 × 10 5 cells/cm 2 ). Cells at low density were incubated in Williams' E medium (Gibco, Invitrogen, Merelbeke, Belgium) supplemented with 10% fetal calf serum (Gibco), 100 units/ml penicillin-100 μg/ml streptomycin (Gibco), 5 μg/ml insulin (Sigma-Aldrich, Bornem, Belgium), 2 mM glutamine (Sigma-Aldrich), and 5 × 10 -5 M hydrocortisone hemisuccinate (Pfizer, Puurs, Belgium). After two weeks of culture, cells were shifted to the same culture medium supplemented with 2% DMSO (Sigma-Aldrich) for two more weeks in order to reach complete differentiation and were maintained for almost two additional weeks. Cells at high density were immediately incubated in the 2% DMSO (Sigma-Aldrich) containing culture medium. Culture media of both low and high density cells were renewed every two or three days [START_REF] Gripon | Infection of a human hepatoma cell line by hepatitis B virus[END_REF][START_REF] Aninat | Expression of cytochromes P450, conjugating enzymes and nuclear receptors in human hepatoma HepaRG cells[END_REF][START_REF] Antherieu | Stable expression, activity, and inducibility of cytochromes P450 in differentiated HepaRG cells[END_REF].

2.2.Isolation of RNA

Total RNA was isolated with the GenElute Mammalian Total RNA Miniprep Kit (Sigma-Aldrich) which directly included a DNase treatment step from HepaRG cells sampled on days 7 and 15 of the proliferation period (n = 6; UNDIFF samples) and on days 1, 3, 10, 12, 15, 16, 18, 19 and 26 of the differentiation period (n = 17; DIFF samples which included cells treated with DMSO to reach complete differentiation and cells at fully differentiated stage, seeded at both densities).

RNA integrity was analyzed with the Agilent 2100 Bioanalyzer™ Automated Gel Electrophoresis System (Agilent Technologies, Diegem, Belgium) according to the manufacturer"s protocol. Its algorithm gives an RNA Integrity Number (RIN) which scales from 1 to 10 with 1 being totally degraded RNA and 10 being totally intact RNA. Sample integrity was determined by the entire electrophoretic trace of the RNA sample [START_REF] Fleige | RNA integrity and the effect on the real-time qRT-PCR performance[END_REF]. This contains the presence or absence of degradation products. As the allocated RIN is independent A c c e p t e d M a n u s c r i p t 6 of sample concentration, instrument and analyst, it is becoming a gold standard for RNA integrity [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF]. The isolated RNA was reverse transcribed with the iScript TM cDNA Synthesis Kit (Bio-Rad, Nazareth, Belgium) followed by cDNA purification with the Genelute PCR clean up kit (Sigma-Aldrich) to remove potential PCR inhibitors [START_REF] De Spiegelaere | Elimination of amplification artifacts in real-time reverse transcription PCR using laser capture microdissected samples[END_REF]. cDNA was finally quantified using a Nanodrop® microvolume spectrophotometer (Thermo Scientific Nanodrop Products, Wilmington, Delaware).

2.3.RT-qPCR analysis

Based on literature search [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Pfaffl | Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations[END_REF][START_REF] Goossens | Selection of reference genes for quantitative real-time PCR in bovine preimplantation embryos[END_REF][START_REF] Cicinnati | Validation of putative reference genes for gene expression studies in human hepatocellular carcinoma using real-time quantitative RT-PCR[END_REF], 12 putative reference genes in the present study were selected for evaluation of their expression profile (Table 1).

Quantative-PCR was performed using either TaqMan hydrolysis probes for TBP, RLP13, ACTB, SDHA, UBC, B2M (PrimeTime qPCR Assay 6-FAM/ZEN/IBFQ, synthesized by Integrated DNA Technology, IDT, Iowa, USA), S18, GAPDH and HPRT1 (TaqMan Gene Expression Assays FAM™ including dye-labeled TaqMan® MGB probe, synthesized by Applied Biosystems, Halle, Belgium) or Perfecta™ SYBR Green FastMix (Quanta Biosciences Inc., Gaithersburg, USA) for YHWAZ, TOP2B and HMBS (synthesized by IDT).

Sequences of primers and probes are outlined in Table 2. All samples were run in duplicate and each run included two negative controls and a serial dilution of a pooled cDNA mix from all samples to calculate the standard curve. The PCR reaction mix consisted of (i) 100 nM of each primer, 12.5 μL Perfecta™ SYBR Green FastMix (Quanta Biosciences Inc.) and 5 µl of DNA in a 25 μl volume adjusted with DNase/RNase-free HLPC water for YHWAZ, TOP2B and HMBS, (ii) 900 nM of each primer, 250nM of probe (both present in 1.25 µl Assay Demand, Applied Biosystems), 12.5 μl TaqMan Universal Master Mix (Applied Biosystems) and 1 or 2 μl of DNA (depending on the tested gene) in a 25 μl volume adjusted with DNase/RNase-free HLPC water for S18, HPRT1 and GAPDH or (iii) 500 nM of each primer, Master Mix (Applied Biosystems) and 1 to 3 μl of cDNA (depending on the tested gene) in a 20 μl volume adjusted with DNase/RNase-free HLPC water for TBP, RLP13, UBC, B2M, SDHA and ACTB. Reactions were performed on the iQ TM 5 Bio-rad system (Bio-rad), with cycling conditions as depicted in Table 2. To verify the accuracy of the YHWAZ, TOP2B and HMBS amplicons, a melting curve analysis was performed after amplification (for YWHAZ en TOP2B from 60° to 94° and for HMBS from 59°C to 94°C, with an interval of 0.5°C and a holding time of 30s).

When possible, probes and primers were designed so that they spanned an exon-exon junction to avoid amplification of genomic DNA.

2.4.Software programs used for statistical analysis and qPCR data processing

Three different algorithms and software programs that calculate the most stable reference genes which are freely available on the Internet were used for statistical analysis and qPCR data processing.

geNorm is a Visual Basic Application for Microsoft Excel which defines the expression stability of a possible reference gene by allocating a so-called M value to each gene in a pool of candidate reference genes. Stepwise exclusion of the reference gene with the least stable expression finally assigns the two most stable genes [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF].

Normfinder is also a Visual Basic applet which assigns a stability value to the candidate reference genes. This robust algorithm, unlike geNorm, adopts a model-based approach to give a score for the two most stable reference genes with the least intra-and inter-group variation. Stability is expressed as a stability value in arbitrary units. Additionally,

Normfinder possesses the ability to discriminate between sample variability and bias between

A c c e p t e d M a n u s c r i p t 8 several groups. The application indeed estimates the intra-and intergroup variance which then describes the stability of the gene expression between the groups resulting in an optimum pair of reference genes [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF].

BestKeeper, another Excel-based tool, determines the variability in expression of a set of reference genes by analyzing quantification cycle (Cq) values and classifying variability by the coefficient of variance (CV) and the standard deviation (SD). To define the most stable reference gene, the software generates a BestKeeper index which finally is compared to each candidate gene resulting in a value for the Pearson correlation coefficient (r) and probability (p) which are then allocated to each candidate reference gene [START_REF] Pfaffl | Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations[END_REF].

Results

3.1.RNA and DNA quantity and quality assurance

The assessed RIN for all 23 samples containing 478 (± 355) ng/µl total RNA fell within a range from 8.6 to 10 indicating the presence of (almost) totally intact RNA, and was hence useful for further qPCR analysis.

Nanodrop®-quantified final cDNA concentrations in all 23 samples varied from 17.8 to 31.9 ng/µl (average = 22.9 ng/µl).

3.2.Statistical and RT-qPCR analysis

PCR-efficiency was calculated from the slope of the standard curves for each gene and fell within the required range of 90-110%. Accuracy of the YHWAZ, TOP2B and HMBS amplicons was proven by a melting curve analysis since one single amplicon was generated. Stepwise exclusion of the reference gene with the least stable expression assigned TBP and SDHA as most stably expressed reference genes in all analyzed samples (DIFF and UNDIFF) (Table 3) as well as in the group of only DIFF samples. In contrast, UBC and YHWAZ showed to be the most stably expressed reference genes in the UNDIFF samples. TBP and SDHA were in the latter group of samples clearly stably expressed as well (Figure 2a-b-c). Finally, by using the pairwise variation between two sequential normalization factors, the software indicated that reliable data for HepaRG can be obtained in RT-qPCR when normalization is performed using a set of the two most stably expressed genes (V2/3 < 0.15) (all samples: TBP and SDHA; UNDIFF samples: UBC and YHWAZ) or the three most stably expressed genes (V3/4 < 0.15) (DIFF samples: TBP, SDHA and UBC). Next to TOP2B and ACTB, the traditionally used reference gene S18, was unstably expressed in HepaRG cells (Table 3).

The solid position of the reference genes assessed by geNorm was appraised by comparison with two other Excel-based applications i.e. Normfinder [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF] and BestKeeper [START_REF] Pfaffl | Determination of stable housekeeping genes, differentially regulated target genes and sample integrity: BestKeeper--Excel-based tool using pair-wise correlations[END_REF]. Normfinder Analysis of our data showed that UBC and TBP were the most stable reference gene in DIFF and UNDIFF HepaRG cells, respectively (Table 3). Distribution of the least stably expressed genes within the DIFF and UNDIFF group of HepaRG samples was slightly different as RLP13, GAPDH and ACTB were the least stably expressed genes in the UNDIFF samples, A c c e p t e d M a n u s c r i p t 10 while TOP2B, B2M and S18 displayed the least stable expression in the DIFF samples (Table 4). When taking into account all analyzed samples (DIFF and UNDIFF), S18, TOP2B and ACTB were again the most unstable genes (Table 3). Following estimation of the intra-and intergroup variance in all HepaRG samples in the present study, the most optimum pair of genes turned out to be TBP and HPRT1.

BestKeeper

Since all genes with a SD >1 are considered inappropriate as reference genes, we could exclude ACTB and GAPDH for analysis when DIFF and UNDIFF HepaRG cells were analyzed both jointly or separately. In addition, when all 23 samples were studied, an SD value of 1.05 was defined for the TOP2B gene and could therefore not be considered as a suitable internal control. Moreover, TOP2B and S18 showed an SD value higher than 1 in the group containing the DIFF samples and were hence not fitting as reference genes in differentiating and differentiated HepaRG cells.

Again, TBP, SDHA and UBC were found to be the most stable reference genes in all 23 samples (Table 3). The BestKeeper analysis also assigned TBP as best reference gene in undifferentiated HepaRG cells, but ranked this gene as fourth stably expressed gene in differentiated cells. However, SD (0,89), r (0,961) and p (0.001) values of TBP were situated closely to the best ranked reference gene RLP13 (SD = 0,76; r = 0.968; p = 0.001) considering TBP also as valuable reference gene when only the DIFF samples are evaluated in the BestKeeper software (Table 4).

When taking geNorm, BestKeeper and Normfinder-analyzed data of all HepaRG samples together, positions for top and bottom ranked genes were fairly constant between the geNorm and BestKeeper software programs. The top three ranked genes included TBP, SDHA and UBC, while ACTB, S18 are TOP2B were listed (nearly) at the bottom of the ranking in both A c c e p t e d M a n u s c r i p t softwares. Rankings differed when using Normfinder, especially for the midfield positioned genes, but this applet also assigned TBP, UBC and SDHA as reliable genes when all 23 samples were analysed.

Discussion

The establishment of a clear-cut set of reliable reference genes is an essential prerequisite for analyzing gene expression using RT-qPCR as shown by several studies [START_REF] Bustin | The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments[END_REF][START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF][START_REF] Ferguson | Impact of Reference Gene Selection for Target Gene Normalization on Experimental Outcome Using Real-Time qRT-PCR in Adipocytes[END_REF][START_REF] Nicot | Housekeeping gene selection for real-time RT-PCR normalization in potato during biotic and abiotic stress[END_REF][START_REF] Erickson | Assessment of normalization strategies for quantitative RT-PCR using microdissected tissue samples[END_REF].

Often "traditional" reference genes, such as S18, are used as internal controls without extensively verifying their validity. As such, the S18 rRNA gene has been used in most studies in which the HepaRG cell line is involved [START_REF] Josse | Long-term functional stability of human HepaRG hepatocytes and use for chronic toxicity and genotoxicity studies[END_REF][START_REF] Vee | Functional expression of sinusoidal and canalicular hepatic drug transporters in the differentiated human hepatoma HepaRG cell line[END_REF][START_REF] Aninat | Expression of cytochromes P450, conjugating enzymes and nuclear receptors in human hepatoma HepaRG cells[END_REF][START_REF] Antherieu | Stable expression, activity, and inducibility of cytochromes P450 in differentiated HepaRG cells[END_REF][START_REF] Cerec | Transdifferentiation of hepatocyte-like cells from the human hepatoma HepaRG cell line through bipotent progenitor[END_REF][START_REF] Dumont | Preferential induction of the AhR gene battery in HepaRG cells after a single or repeated exposure to heterocyclic aromatic amines[END_REF][START_REF] Lambert | Reproducible chemical-induced changes in gene expression profiles in human hepatoma HepaRG cells under various experimental conditions[END_REF]. Next to human primary hepatocytes, this hepatocyte model seems to be a very suitable tool for in vitro toxicity testing of xenobiotics and even has some additional advantages [START_REF] Guguen-Guillouzo | Stem cell-derived hepatocytes and their use in toxicology[END_REF][START_REF] Guillouzo | Evolving concepts in liver tissue modeling and implications for in vitro toxicology[END_REF]. Indeed, the cell line can be used either undifferentiated (progenitor cells) or differentiated (hepatocyte-like and biliary-like cells) making it possible to distinguish between compounds that are directly toxic and those that need activation to be hazardous to either cell type. It generally acquires both the metabolic capacity of freshly isolated human hepatocytes as well as the potential to grow indefinitely like cancer cells [START_REF] Guguen-Guillouzo | Stem cell-derived hepatocytes and their use in toxicology[END_REF][START_REF] Guillouzo | Evolving concepts in liver tissue modeling and implications for in vitro toxicology[END_REF]. This commercially available cell line is indeed increasingly adopted at the academic as well as industrial level. Hence, when performing RT-qPCR experiments, it is crucial to include appropriate normalization controls due to their progressively apparent significance as the amount of biological studies applying this methodology has increased.

To diminish the menace of erroneous results because of the use of invalid reference genes, we As mentioned by Vandesompele et al. [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF], mRNA expression data are much more reliable and accurate when they are normalized using the geometric mean of multiple reference genes.

As such, geNorm analysis demonstrated that, for the current study design, two valuable reference genes (TBP and SDHA) are required for normalization when analyzing expression levels in HepaRG cells at proliferating and differentiating or differentiated stage (i.e. all samples and UNDIFF samples). However, when mRNA expression was analyzed in merely the group of differentiating and differentiated cells (DIFF samples), geNorm assessment showed that adding the 3 rd candidate reference gene (TBP next to UBC and YHWAZ) to the normalization factor, had the largest influence on diminishing variability. This undoubtedly underscores the previously stated importance of selection of a sufficient set of reference genes because of the variability within a single cell type across different experimental treatments.

GeNorm additionally showed that S18, commonly used as single internal control in HepaRG studies, was nearly the least stably expressed reference gene in HepaRG cells at all stages.

This finding was confirmed by Normfinder and BestKeeper analyses. Collectively, these data show that S18 is not suitable for normalization in the HepaRG cell system -and definitely not as single reference gene -resulting in possible misinterpretation of qPCR data.

In addition, current analyses showed that geNorm and BestKeeper both ranked the examined candidate genes generally in comparable order, certainly those positioned in the top and bottom ranking. Normfinder, however, provided a slightly different list of ranked genes. It defined for example HPRT1 as a reference gene of the most optimum pair next to TBP, whereas the HPRT1 gene was only positioned as the 7 th and 8 th most stable reference gene in the geNorm and BestKeeper analyses, respectively. This may be explained by the fact that Normfinder evaluates the candidate as to whether it shows systemic variation across the sample subgroups and hence does not only take overall expression variation into account [START_REF] Andersen | Normalization of real-time quantitative reverse transcription-PCR data: a model-based variance estimation approach to identify genes suited for normalization, applied to bladder and colon cancer data sets[END_REF].

Nonetheless, UBC and TBP were also highly ranked as potential reference genes by this applet, while TOP2B and S18 were considered less valuable or even useless as reference genes in the HepaRG cell line.

Discrepancy was also seen for some genes when DIFF and UNDIFF samples were separately analysed, for example for the B2M gene and the ACTB gene, which were differently ranked within both groups by the three software programs and Normfinder, respectively. This finding again emphasizes the need for testing several candidate reference genes in HepaRG cells since some reference genes are differently regulated across multiple experimental conditions. This may be, for example, the case for studies wherein HepaRG cells are treated with different inducers of biotransformation phase I and phase II enzymes.

Although as few as two or three reference genes can serve as a good standard [START_REF] Erickson | Assessment of normalization strategies for quantitative RT-PCR using microdissected tissue samples[END_REF][START_REF] Derveaux | How to do successful gene expression analysis using real-time PCR[END_REF] which was in our study demonstrated by geNorm analysis too, we would suggest, when setting up RT-qPCR experiments in HepaRG, to use an initial set of minimum six to eight candidate reference genes which have been confirmed in the current study to be the most stably expressed. When analyses of the three software programs are combined, this set of genes contains TBP, UBC, SDHA, RLP13, YHWAZ, HMBS, B2M and HPRT1. From these candidates, minimum two reference genes can then be selected for normalization in a specific experimental design with HepaRG cells.

In conclusion, while TBP was proven to be overall the most stable reference gene, TOP2B, ACTB and S18 were considered less valuable as internal control in HepaRG cells. In general, we propose to test a set of at least six to eight reference genes (TBP, UBC, SDHA, RLP13,

Table

A c c e p t e d M a n u s c r i p t

A c c e p t e d M a n u s c r i p t 7 250

 7 nM of probe (both present in 1 µl 20X Prime Time Assay, IDT), 10 μl TaqMan Universal

Figure 1

 1 Figure1shows the cycle threshold value (Cq-value) of the 12 tested reference genes which

M a n u s c r i p t 12 (

 12 validated 12 potential reference genes, using geNorm, BestKeeper and Normfinder algorithms, suitable for RT-qPCR profiling experiments in this promising in vitro liver-based model during proliferation (UNDIFF), during differentiation and at fully differentiated stage A c c e p t e d DIFF). Cultivation was performed exactly as described in the literature [2, 9, 10].

A c c e p t e d M a n u s c r i p t 13

 13

A c c e p t e d M a n u s c r i p t 14 YHWAZ

 14 , HMBS, B2M and HPRT1) of which then the most stable genes can be picked out for further normalization in future HepaRG expression studies. As such, reproducibility and validity of RT-qPCR experiments in this in vitro liver-based model may notably improve.

Figure 1 .

 1 Figure 1. Cq-value range of the reference genes.

Figure 2 .

 2 Figure 2. Reference gene mRNA expression stability according to geNorm (a) for all analyzed

Figure A c c e p t e d M a n u s c r i p t

  Figure

Figure 1 .

 1 Figure 1.

  Figure

Figure 2b.

  Figure 2b.

Figure A c c e p t e d M a n u s c r i p t

  Figure

Figure 2c.

  Figure 2c.

Table 1 .

 1 Twelve selected candidate genes, their encoded proteins and biological function.YHWAZ 14-3-3 protein zeta/delta Protein involved in cell signaling, regulation of cell cycle progression, cytoskeletal structure, and transcription [35] TOP2B topoisomerase (DNA) II b Enzyme controlling and altering the topologic states of DNA during transcription [36]HMBSHydroxymethylbilane synthase Heme synthesis, porphyrin metabolism[START_REF] Derveaux | How to do successful gene expression analysis using real-time PCR[END_REF]

	Abbreviated gene name	Encoded protein

Table 3 .

 3 Ranking of candidate reference genes in all 23 analyzed samples based on geNorm, Normfinder and BestKeeper analyses.

	Gene		geNorm		BestKeeper		Normfinder
		M value	Ranking	r	p	Ranking	Stability value	Ranking
	GADPH	0.896	8	Nd	Nd	Nd	0,888	9
	HPRT1	0.854	7	0.864	0.001	8	0,300	2
	S18	1.381	11	0.694	0.001	9	0,915	10
	ACTB	1.049	10	Nd	Nd	Nd	1,213	12
	B2M	0,903	9	0.906	0.001	6	0,499	4
	SDHA	0,730	2	0.961	0.001	2	0,534	5
	UBC	0,745	3	0.953	0.001	3	0,397	3
	RLP13	0,746	4	0.936	0.001	4	0,876	8
	TBP	0,712	1	0.986	0.001	1	0,159	1
	YHWA Z	0,756	5	0.900	0.001	7	0,551	6
	TOP2B	1,459	12	Nd	Nd	Nd	1,134	11
	HMBS	0,763	6	0.910	0.001	5	0,717	7

M value, expression stability measured in geNorm (should be < 1.5); Nd, Not determined (excluded for final BestKeeper analysis because of SD >1); r, correlation coefficient (BestKeeper) (should be close to 1); p, probability value (BestKeeper); Stability value (Normfinder) should be as low as possible.

Acknowledgements

This work was supported by grants from the Industrial Research Funds, (IOF), Vrije Universiteit Brussel, Belgium. Special thanks go to the European Union (FP6 project carcinoGENOMICS). We also thank Prof. P. Simoens (Ghent University) for his critical comments.

The authors thank you for reviewing the manuscript entitled "Critical selection of reliable reference genes for gene expression study in the HepaRG cell line", by Ceelen L. et al.

We believe that it is an excellent proposal to check whether conventional inducers of phase I and II enzymes would influence the expression profile of the 12 candidate reference genes. However, the main message of this manuscript, as stated in the abstract and discussion, was to propose a set of approximately eight reference genes which should be validated in every new experimental design using the HepaRG cell line. We also wanted to investigate whether the commonly used S18 gene was a valid reference gene, which turned out not to be. We believe that it is of utmost importance to analyze this set or reference genes when, as you correctly suggested, HepaRG cells are subjected to various experimental conditions. Following validation of this set of genes, the most stable ones can then be selected under these specific conditions. So, the take home message of the paper is that researchers should not use S18, or any other single reference gene for normalization in HepaRG gene profiling studies as such, but first to validate a set of genes and then pick out the most stable ones (a minimum of two reference genes is required as indicated by the geNorm applet [START_REF] Vandesompele | Accurate normalization of real-time quantitative RT-PCR data by geometric averaging of multiple internal control genes[END_REF]; MIQE guidlines).

In addition, in the current set, TBP was indeed the most stable reference gene, but the authors clearly declared not to use TBP as single reference gene (but as stated before to use a minimum of two genes). Depending on the experimental design (for example HepaRG cells treated with different inducers of Phase I and II enzymes), TBP may be again the most stable reference gene, but it should be validated against the other candidate reference genes in every new experimental setting, as you correctly observed.

The authors want to mention in addition, that at the current moment, experiments using HepaRG cells under different experimental conditions are running wherein the proposed set of candidate genes is being validated as well. These experiments have been set up in the framework of a large project and results will be published in the future.

The authors regret that the message of the paper was not entirely clear, but we hope by adapting the final conclusion in the abstract and by adding your proposal to the discussion, that the final point of the current paper is now clear. We also want to thank the reviewer for his/her excellent proposal. To additionally underscore the take home message to the reader, we have changed the last sentence in the abstract (i.e. the main conclusion) and added a sentence to the discussion for further clarification. Moreover, due to impracticalities, it would not be possible to perform additional experiments as suggested. We also have added a letter to the reviewer in order to clarify the ultimate message of this paper.

The authors did a final review of the paper including a spelling and lay-out check (see below). We also want to confirm that all suppliers and their location are provided in the manuscript text for every chemical used in this study. Finally, we noticed that the melting curve analysis (for YHWAZ, TOP2B and HMBS) was not mentioned in the original text, although this analysis was performed. Therefore, to complete the RT-qPCR analysis, this information has been added to the revised text.

Textual changes include:

1. Page 3: These cells have been used for long-term repeated exposure to evaluate the potential toxic effects of different chemical entities [START_REF] Guguen-Guillouzo | Stem cell-derived hepatocytes and their use in toxicology[END_REF][START_REF] Josse | Long-term functional stability of human HepaRG hepatocytes and use for chronic toxicity and genotoxicity studies[END_REF] and. -> "and" has been removed. 2. Page 5: 100 units/ml penicillin,100 μg/ml streptomycin (Gibco)-> 100 units/ml penicillin-100 μg/ml streptomycin (Gibco) 3. Page 6 and 8: Nanodrop has been replaced by Nanodrop® 4. Page 7: iQ5 Bio-rad system has been replaced by iQ TM