Roger Collet
email: rcollet@inrets.fr.

Matthieu De Lapparent
email: lapparent@inrets.fr.

Laurent Hivert
email: hivert@inrets.fr.

Addiction to car use and dynamic elasticity measures in France

Keywords: Transportation, Car use, Consumption, Addiction, Panel, GMM. JEL Classification: C23, D12

 is applied on a panel dataset, drawn from the French "Car Fleet" survey over the period 1999-2001. Importantly, the estimates show that the assumption of addiction to car use cannot be rejected. Furthermore, the model yields realistic kilometric-price and income elasticities of household automobility, for both the short and the long runs.

Introduction

Contribution of mankind activity to climate change and the related ecological issues are subject to tempestuous political and scientific discussions. Whoever may be right about the exact assessment of causes and consequences, French policy makers have decided to set more and more ambitious environmental targets, and measures are implemented to converge towards a more sustainable development. Voted in 2005 by the French Parliament, the energy law aims at dividing by 4 the 1990 level of greenhouse gases at the horizon 2050. In a more global perspective, we refer the reader to the debated results since the 1979 world climate conference in Geneva up to the 2009 climate change conference in Copenhagen. Nowadays, there is a large international consensus about the emergency to make the economic growth cleaner. However, it is noteworthy that there is no general agreement about the way to manage it. Transportation is a sector for which a particular effort is deemed essential. Industrial strategies and policies concern both supply and demand sides. As it regards the latter, the subtopic of travel demand focuses on behavioural changes and incentives to turn towards more sustainable means of travelling. A special attention is paid to emission of greenhouse and other toxic gases but not only. Indeed, the perspective of global peak oil and oil scarcity pleads for the necessity to find new and/or more efficient ways to travel, as it should entail an increase of fuel prices in a near future. Therefore, the question of car use intensity emerges, given that current technology is mainly based on fossil resources. In a context of highly volatile and increasing fuel prices, it is of great importance to understand and to quantify to what extent households adapt their behaviours in terms of car ownership and use. Although not that much new, this remains a major topic of research. Car provides users with a larger control of space and time, and allows a better access to jobs, leisure places, health care, public amenities... On the other hand, car users may develop a kind of dependence, as pointed out by Dupuy (1999, p.1):

"… automobile dependence means that as individuals, we cannot live without cars, just as a smoker cannot live without cigarettes, and a drug addict without drugs."

Below, Wickham (2002, p.16) also compare car use with an addictive consumption of drugs. Our microeconomic reading is proposed in parentheses.

"Car (use) dependence can be understood through the metaphor of drug dependency (addiction). Heroin or even nicotine addiction is in part a matter of (rational) choice. I choose to shoot up, I choose to smoke a cigarette. But as I continue to do this, my body (utility function) changes, it becomes restructured, it needs the drug (addictive good), it cannot do without it. Furthermore, the need (addictive good marginal utility) escalates the body (optimal bundle of goods) requires more and more of the drug (addictive goods)... the same applies to car (use) dependency."

In the present article, we propose a microeconometric analysis of the annual mileage travelled by French households with their personal cars. In the literature, this mileage is sometimes referred to as household "automobility", and we will use indifferently both expressions all along the paper. A two-step structure is applied: after considering the dichotomous choice of households to own cars or not, their automobility is modelled. An objective focuses on measuring price and income elasticities of car use. Moreover, the panel layout of the data, based on three annual waves (1999)(2000)(2001) of the French "Car Fleet" survey, enables to use dynamic specifications and to derive short and long-run effects. To feature the dependence of households to car use, the myopic and rational addiction models [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] are investigated. By applying these models, generally used to describe the consumption of cigarettes, alcohol or drugs, the previous authors' assertions can be tested from the microeconomic point of view. To our knowledge, this article provides the first application of the rational addiction model to describe car use behaviour, while "automobile dependence" is a major topic of research in transportation. The rest of the article is fashioned as follows. Section 2 presents the theoretical framework: the rational addiction model and its properties are examined. Section 3 discusses the econometric specification and the estimation method. Section 4 provides a description of the 1999-2001 French "Car Fleet" panel dataset. Section 5 reports and debates the estimation results. Conclusions are drawn in a last section and further research tracks are also suggested.

Microeconomic framework

2.1

The rational addiction model

Description

Since [START_REF] Becker | A theory of rational addiction[END_REF], a consumer is said to be addict to a good if, all the same, an increase in his past consumption yields a significant rise in his current consumption. [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] developed an addiction model in which the individual current utility level t U depends on the consumed quantity of two goods: a quantity t X of a composite good X , and a quantity t C of an addictive good C . The current utility also depends on a set of individual characteristics t e related to lifecycle and potentially unobserved. The addictive and composite goods are different in that t U also depends on the past consumed quantities of C . Following the authors, these quantities are accumulated into an addictive capital stock t S , which is supplied at each current period by the past level of consumption 1 t C -. Thus, the current utility of a consumer also depends on this stock. The most commonly used expression to define it simply writes:

1 t t S C - = .
The consumer looks to maximize his discounted utility over an infinite lifetime horizon:

( )

1 1 1 , , , t t t t t t t Max B U C C X e ∞ - - = ∑ , (1)
where ( )

1 1 B ρ - = +
is the discount factor, ρ being the inter-temporal rate of substitution. Besides, it is assumed in the model that the composite X is the money and the interest rate of the economy is equal to ρ . Maximization of the consumer's utility is subject to an initial condition regarding the addictive good, and it is constrained by the inter-temporal budget equilibrium:

0 0 C C = ;
( )

1 0 1 t t t t t A B X PC ∞ - = = + ∑ , ( 2
)
where 0 A is the net present value of the consumer's wealth, t P is the nominal price of the addictive good at date t . Let the utility function be concave and quadratic in the arguments ( )

2 2 1 1 1 2 2 1 , , , 2 2
= + + + , (4)
where:

( ) ( ) ( ) (
)

1 2 CC XX CX SS XX XX CS B θ α α α α α α α - = - - + . ( 5
)
For CS α strictly positive, the past and current consumptions of C are said to be complementary. In this case, the current marginal utility of the addictive good, say ' t C U , is an increasing function of

= = + + + + . ( 6
)
The higher the level of 1 t Cand the value of CS α , the higher the marginal utility of the addictive good. By analogy with "learning by doing", the consumer enjoys all the more the consumption of the addictive good as he "practiced" it in the past, and as the "learning speed" ( CS α ) is high.

The temporal complementarity of the consumptions of C is the origin of addiction. It implies 0 θ > in (5). The larger the estimated value of θ , the greater the level of addiction. The static and autoregressive consumption models are particular cases of the dynamic demand function (4). For 0 θ = , the demand function does neither depend on the past nor on the future levels of consumption and the static case emerges. When 0 θ > and 1 / (1 ) B ρ = + drops to 0 (that is, for an infinite preference for the present), the demand function (4) takes the form of a firstorder autoregressive process. There, the individual is not a forward looking agent. [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] defines such a behaviour as myopic addiction. In any other situation where θ and B are strictly positives, the addiction behaviour is said to be rational. As both these parameters can be estimated, it is possible to test the most relevant formulation from the empirical perspective.

Testing for rational addiction

Testing the rational addiction theory has been performed on various data that pertain to different topics of research. Many applications dealt with drug consumption to tackle public health issues. For instance, [START_REF] Baltagi | The Econometrics of Rational Addiction: The Case of Cigarettes[END_REF], [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF], [START_REF] Gardes | Evidence on Addiction Effects from Households Expenditure Surveys: the Case of the Polish Panel[END_REF], and Tiezzi (2004) used addiction models to explain the consumption of tobacco. Grossman et al. (1998), [START_REF] Bentzen | Rational Addiction and Alcohol Consumption: Evidence from the Nordic Countries[END_REF], [START_REF] Baltagi | Rational addiction to alcohol: panel data analysis of liquor consumption[END_REF], and [START_REF] Lalla | Rational addiction for alcohol: a panel study of Norwegian beer drinkers[END_REF] applied them to model the demand for alcohol. [START_REF] Van Ours | The Price Elasticity of Hard Drugs: The Case of opium in the Dutch East Indies, 1923 -1938[END_REF], and Grossman and Chaloupka (1998) tested the relevance of addiction models to explain the consumption of opium and cocaine.

As mentioned by [START_REF] Becker | A theory of rational addiction[END_REF], there is no need to express a biological dependence to be considered as addicted to a good. There are other topics of application than drug consumption. For instance, [START_REF] Mobilia | Gambling as a Rational Addiction[END_REF] tested for addiction to gambling, [START_REF] Cawley | Rational Addiction, the Consumption of Calories, and Body Weight[END_REF] focused on the consumption of calories, [START_REF] Villani | Rational Addiction in the Arts[END_REF] dealt with addiction to art, [START_REF] Cameron | Rational Addiction and the Demand for Cinema[END_REF][START_REF] Sisto | Rationally Addicted to Cinema and Tv? An empirical investigation of Italian consumers[END_REF][START_REF] Sisto | Rationally Addicted to Cinema and Tv? An empirical investigation of Italian consumers[END_REF] applied addiction models to the demand for cinema. To our knowledge, the rational addiction model has not been applied yet to describe car use, while "automobile dependence" is a major topic of transportation research.

One needs however to be careful about the results of such models. Actually, whereas the effects of past and future consumptions on the current one are often found to be significant, the resulting inter-temporal rates of substitution are sometimes not convincing. [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] reported some rates ranging from 56 up to 223 per cent. Tiezzi (2004) also found unlikely values. [START_REF] Baltagi | Rational addiction to alcohol: panel data analysis of liquor consumption[END_REF], Grossman et al. (1998), even reported some negative rates. All in all, a rule of thumb would be to accept the assumption of rational addiction on the basis of a "reasonable" rate1 in addition to the statistical significance of past and future consumptions.

Measures in the rational addiction model

Although the core specification to test for the addiction hypothesis remains the model ( 4), it is often generalized to incorporate additional variables 2 . In the present approach, we consider that the current covariates are of two types: economic-related variables it E (nominal price and income for example) and else-related variables it S (socio-demographic and geographic variables for example). From the modeller's point of view, it E and it S may both include observed and unobserved factors. These latter are all summarized by an error term it ε . Let our specification be:

1 1 0 1 it it it it it it C C BC S E θ θ α α ε - + = + + + + , (7)
where the subscripts i and t identify the individual and the period. The impacts on the current consumption resulting from a variation in the past and future consumptions can be deduced from the characteristic roots of the homogeneous equation of ( 7). They write:

1

1 1 4 ² 2 B θ ϕ θ -- = , 2 1 1 4 ² 2 B θ ϕ θ + - = . (8)
In ( 8), 1

ϕ and (Becker, 1996, p.113):

( ) ( ) 1 / 1 2 1 CT k k C Ek E e C α θ ϕ ϕ = × - , (9) ( )( ) (
)

1 / 1 2 1 1 LT k k C Ek E e C α θ ϕ ϕ - = × - - . (10)
Logically, the elasticities that stem from the myopic model are particular cases of ( 9) and ( 10) for 1 0 ϕ = and 1 2 ϕ θ -= .

Modelling and estimation difficulties

Notations

Unless explicitly stated, the following notations will be used all throughout the rest of the paper. In addition to the subscripts i and t which identify respectively the household and the period, we consider also a subscript v which is related to a specific car owned by i at date t.

These three subscripts are used all together to index the following variables:

• KM: the annual mileage converted into kilometres;

• FE: the average fuel efficiency in litres per 100 kilometres;

• DP: the price per litre of diesel oil in Euros;

• PP: the price per litre of premium-petrol in Euros;

• IP: the average of the two latter prices when information about the type of fuel used by the engine is not known;

• KMP: the kilometric price, or more precisely, the fuel operating cost per 100 kilometres. It is defined as

ivt ivt ivt KMP FE FP = ×
, FP being the price of the type of fuel (DP, IP, or IP);

• AD: the age of the car if it is diesel-powered;

• AP: the age of the car if it is petrol-powered;

• AI: the age of the car if the type of fuel used is not known;

• AN: a dummy variable that states whether the age of the car is not known.

Additional variables that describe the characteristics of the household will also be used. They are detailed in the next subsection.

3.2

The econometric model

The annual mileage of a car v that is owned by household i at date t is modelled as:

0 1 2 3 4 5 ( ) , ivt s it ivt ivt s ivt ivt ivt ivt KM R s KMP AD AP AI AN β β β β β β ε = = + + + + + + ∑ 1 (11)
The error term ivt ε is assumed to be drawn from a zero-mean normal distribution. The model intercept 0

β is here differentiated according to the household residential location: either Paris ). The automobility of a household i at date t is got by summing the mileages of the cars it owns at this date. Thus aggregated, the following model emerges:

0 1 2 3 4 5 ( ) , it ivt s it it ivt ivt v s v v ivt ivt ivt ivt v v v v KM KM R s NC KMP AD AP AI AN β β β β β β ε = = = + + + + + + ∑ ∑ ∑ ∑ ∑ ∑ ∑ ∑ 1 (12)
where NC refers to the number of cars.

The previous specification is enlarged by introducing additional variables to allow for a better control of households heterogeneity: the household annual income in Euros, three dummy variables describing the age class of the household head ([18-39], [40-65], >65), the number of adults (except the head), the number of working adults, the number of women, the number of driving-license owners, and the number of children in the household. Finally, according to the type of addiction model (myopic or rational) to be estimated, the specification also includes the past and future household automobility ( 1 t KM ± ) as explanatory variables.

Selectivity, heteroskedasticity and endogeneity

The addiction models require that the dependant variable is not censored. The fact is that the survey from which we have drawn the data contains households that did not own cars. The annual mileage for these non-motorized households is zero, corresponding in microeconomics to a corner solution. The automobility models have been estimated using the subsample of households that declared to own at least one car in 2000. However, excluding the nonmotorized households may result in a sample selection problem.

We have accounted for this likely selection bias by applying a two-step estimation procedure [START_REF] Heckman | Sample Selection Bias as a Specification Error[END_REF]. In a first step, we have estimated a dichotomous Probit model that explains household car ownership in 2000. We have then used the results to estimate the inverse Mills ratio λ , say λ . The latter has in turn been used as an independent variable in the household mileage equation, which has been estimated using the subsample of car owner households in 2000. This corrects the potential bias by capturing the correlation between the ownershipbased selection and the automobility models. Testing for selection bias is easily carried out by checking whether the estimated coefficient that weighs λ , say λ β , is statistically different from zero. This correction method of sample selection is not without some difficulties. Actually, the introduction of λ into the automobility model generates heteroskedasticity, as it makes the variance of the related error depend on the covariates used in the ownership Probit model [START_REF] Heckman | Sample Selection Bias as a Specification Error[END_REF]. Besides, the use in a specification of pre-estimated variables ( λ in our case) as determinate covariates generally leads to underestimate the estimator variance3 . In such circumstances, the correction provided in [START_REF] Murphy | Estimation and inference in two-Step econometric models[END_REF] can be applied. Sample selection is not the only reason why heteroskedasticity has to be taken into account.

The structure of the proposed model is also a source. For instance, stating that the error terms ivt ε are identically and independently distributed is not a convincing assumption. Indeed, a negative shock on the error of a household's car (as a breakdown) can induce a positive impact on the error of other household's cars. Typically, this is designing a correlation pattern of errors between the cars of a same household:

' cov( , ) 0 ivt iv t ε ε ≠ . Besides, random shocks,
caused by particular traffic conditions or public transport supply for example, can be unequally passed on to the car mileages of a household. They should also depend on the location and the ownership level:

' ' var( ) var( )

ivt i v t ε ε ≠ for i i′ ≠ or v v′ ≠ .
At last, the summation of car mileages to compute household automobility is also producing heteroskedasticity, related to the car ownership level. Therefore, an estimation method allowing for a complex form of heteroskedasticity must be applied. By involving simultaneously both the lagged and the forwarded dependant variables as covariates, the dynamic specification of the rational addiction model makes them necessarily endogenous, even assuming the temporal independence of individual errors. Moreover, the errors are likely to be serially correlated due to an unobservable time invariant heterogeneity factor i η , assuming that v ivt

Estimation strategy and tests

A solution resides in turning to estimators that use instrumental variables. Among those existing, the 2SLS estimator has almost but not all the desired properties. Although convergent, it is not consistent in the presence of heteroskedastic error terms. It may however be used to test for existence of heteroskedasticity (Breush and Pagan, 1979). If the test rejects the homoskedasticity assumption, the Generalized Method of Moments (GMM) should then be applied. Proposed by [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF], this method generalizes other estimators such as the OLS and the 2SLS estimators. A condition of application of the GMM is however to use a set of "good" instruments. They need to be orthogonal to the estimation residuals and enough correlated with the endogenous variables to be instrumented. Both properties have therefore to be examined. To that extent, tests proposed by [START_REF] Hansen | Large sample properties of generalized method of moments estimators[END_REF] and [START_REF] Bound | Problems with instrumental variables estimation when the correlation between the instruments and the endogeneous explanatory variable is weak[END_REF] are implemented. The readers are referred to [START_REF] Baum | Instrumental variables and GMM: Estimation and testing[END_REF] for a detailed discussion.

Data and descriptive statistics

Data source

Data are drawn from the French "Car Fleet4 " panel survey, which is achieved annually since 1984 by the private pooling institute TNS-Sofres5 . The survey aspires to a better knowledge of several dimensions of the automobile demand, especially car ownership and use. It depicts with a great level of details many of the attributes of the cars owned by the households, as well as many of their characteristics. A nationally representative sample of 10 000 households is surveyed each year. The panel is rotating: about one third of the sample is renewed each year. Such a methodology allows for a longitudinal follow-up of some households for at least 3 years. In the present paper, we focus on households that were surveyed over the period 1999-2001.

Descriptive statistics

We have identified 3010 households continuously present in the waves 1999, 2000 and 2001 of the survey. On annual averages, slightly less than 20 per cent of them have no car, about 50 per cent have one car, slightly more than 25 per cent have two cars and about 5 per cent have three cars or more. The average automobility of households is monotonically decreasing over the considered period, from 15 610 kilometres in 1999 down to 14 826 kilometres in 2001. Excluding households without a car to cancel out the decision related to car ownership, we observe the same decreasing trend: the average household mileage decreases also monotonically from 19 279 kilometres to 18 189 kilometres over the period. Table 1 reports the annual descriptive statistics about the characteristics of the sampled households. In 1999, a total of 3552 cars are reported by the households in the panel. It is a total of 3576 cars in 2000 and a total of 3605 cars in 2001. Based on this sample, Table 2 provides the descriptive statistics regarding car attributes and mileages. In 2000, the average car is almost 6.8 years old. According to the engine type, petrol-powered cars are about two years older on average than diesel-powered cars. This difference results from the very dynamic trend of global dieselization of the French car fleet [START_REF] Hivert | Dieselisation and the 'new dieselists' behaviour: recent developments in the French car fleet[END_REF]. In 1980, diesel cars represented less than 5 per cent of the total fleet in France. This share has continuously increased to reach about 15 per cent in 1990, 30 per cent in 1995, 35 per cent in 1999 and 40 per cent in 2001 6 . In accordance with these figures, the proportion of diesel cars is also increasing in our data, from 35 per cent in 1999 to 38 per cent in 2001. Dieselization explains partially the improvement in time of the average energy efficiency of vehicles, since diesel cars are less fuel-consuming than petrol cars of about 0.9 litres for 100 kilometres (Table 2). It also derives from the improvement of the fuel efficiency for both types of car over time: globally in the data, the average vehicle consumed 7.33 litres of fuel for 100 kilometres in 2001 to 7.44 litres initially in 1999. After a period of low fuel prices during the 1990's, the year 2000 marked an episode of significant raise: in 1999, due to the decision of various oil-producing countries (including OPEC) to limit the production, the price of crude oil barrel has increased. On annual averages, the price per litre of premium-petrol in France raised from EUR 0.98 in 1999 to EUR 1.14 in 2000 (from EUR 0.69 to EUR 0.85 for diesel-oil). But this increase has been short-lived, because during 2000, the production of oil has increased again, resulting in a decline of fuel prices for the following year: in 2001, the price for one litre of premium-petrol dropped to EUR 1.09 (EUR 0.80 for diesel-oil). However, the 1999-2001 whole result still represents an increase in price for both premium-petrol (+11%) and diesel-oil (+16%) (Table 3). These differences in fuel price and energy efficiency partly explain the more intensive use of diesel cars. These ones have covered 17 085 km on average in 2000 to "only" 10 412 km for petrol cars. Between 1999 and 2001, the average mileages have decreased by 1600 km and 840 km for diesel and petrol cars respectively (Table 2). Thus, the decrease has been higher for diesel-powered cars than for petrol-powered cars. A reason is that the diesel-oil price per litre has increased faster than the petrol price.

Results

Table 4 reports the estimates of the selection model, which is the top layer of our modelling approach. It models the household probability to own at least one car in 2000. These estimates are then used to compute the correction factor λ , which is introduced in the automobility model to control for selection. The estimates of the addiction models are reported in Table 5.

Bound (p-value)

F(22, 2590) KMt-1: 3.80 (0.00) F(22, 2590) KMt-1: 3.80 (0.00) KMt+1: 4.51 (0.00) Note: subsample of 2631 households which described at least one car in 2000. Denominations of variables are in section 3.0. Dependant variable: KM in 2000. GMM estimation with 11 excluded instruments for 1999, 11 for 2001, 18 current instruments included in the specification. T-statistics in parentheses. The Breush-Pagan statistic is computed from the residuals of models estimated by 2SLS.

Fit properties

Both myopic and rational addiction models present good fit properties. The ² R statistics are respectively 0.87 and 0.90, and the Fisher statistics show that the set of explanatory variables is relevant to explain the dynamics of household automobility. The Breush-Pagan statistic after estimating the models by 2SLS rejects the homoskedasticity assumption, and justifies resorting to the GMM estimator as IV technique. The instruments that were used to implement the GMM are the set of current exogenous covariates ("included" instruments) and a set of past and future household characteristics ("excluded" instruments). In both models, the Hansen test concludes to acceptance of the null hypothesis of orthogonality between residuals and instruments. Besides, the Bound test accepts the alternative hypothesis of joint significance of the excluded instruments to explain the endogenous covariates.

Selectivity

Both models agree not to accept the significance hypothesis of the correction term λ . Formally, the null hypothesis Ĥ0 : 0 λ β = cannot be rejected at the 5 per cent level, meaning that the selection of car owner households in 2000 to estimate the automobility models has not been a source of bias. Under H0, the estimator variance corrected according to the Murphy-Topel method amounts to the variance without correction 7 . Thus, the uncorrected estimator variance should not differ statistically from the corrected variance. Moreover, whichever the addiction model, the Hausman test leads to accept the null hypothesis of equal estimates between the unconstrained model and the model constrained by ˆ0 λ β = . Therefore, using the unconstrained estimates and the uncorrected estimator variance (from Table 5) for inferential purposes is legitimate.

Addiction and inter-temporal rate of substitution

While the results substantiate the addictive behaviour regarding automobility, the rational addiction model emerges as the most relevant to describe the sort of car dependence of households captured by the data. Indeed, the myopic model confirms the significant effect of the past annual mileage on current automobility. But on the other hand, the rational model rejects behavioural myopia as the parameter that weighs the forwarded annual mileage is also statistically significant. Therefore, households are forward-looking agents in setting their current automobility. Furthermore, the rational addiction model yields an inter-temporal rate of substitution of about 17 per cent, which is a rather plausible value.

Price and income elasticities of automobility

The estimates related to the kilometric price have the expected negative signs in both models, but their levels of significance slightly differ. It is significant at the 1 per cent level in the myopic model while it is just significant at the 5 per cent level in the rational model. Elasticity measures of household automobility to the kilometric price are reported in Table 6.

The myopic and rational addiction models agree to measure a short-run price-elasticity of 7 Referring to the correction formula in Greene (2003, p.510), the corrected estimator variance ( * Relatively, the short-run elasticity obtained from this model (-0.22) might seem high, meaning a fast rate of convergence to the long-term equilibrium. In other words, French households adapt quickly their automobility to a change in the kilometric price. This conclusion can also be drawn from Graham and Glaister (2002, p.22, fig.1). Indeed, these authors reported the price-elasticities of the demand for gasoline for a set of occidental countries: comparatively, France has one of the highest sensitivities in the short run and one of the lowest for the long run. The short-long ratio is about 3/4, while it is clearly below 0.5 for the other comparable countries (Germany, United Kingdom, Austria and Canada). Table 8 reports the sensitivity of car annual mileage to fuel price by type of car. The elasticity of mileage for petrol-powered cars with respect to premium-petrol price is estimated at -0.32 in the short run and at -0.52 in the long run. Regarding diesel-powered cars, these elasticities with respect to diesel-oil price are respectively -0.13 and -0.21 8 . Thus, the use of petrol cars is about 2.5 times more sensitive to fuel price variations than diesel cars, both in the short and the long runs. The estimated coefficient that pertains to households' annual income is significant and has the expected positive sign in both models. Table 6 reports also the short and long-run elasticities of household automobility with respect to income. Whatever is the considered addiction model, these elasticities take similar values: about +0.10 in the short run and about +0.16 in the long run. Thus, an increase in household income entails a proportionally lower increase in automobility: in the economic sense, the kilometer proves to be a normal good. Using the results from the rational addiction model, a permanent increase by EUR 1000 in household 8 Example: in 2000, the average fuel consumption for a diesel car was 6.80 litres of diesel-oil for 100 km, the average mileage for this type of car was 17 085 km (Table 2), and the price for one liter of diesel-oil was EUR 0.8461 (Table 3). Moreover, the long-run marginal effect of the kilometric price on automobility is estimated at -623 km by car (Table 7). For diesel cars, the long-run elasticity of mileage to the diesel-oil price is valuated at: (-623×6.80)×(0.8461/17085)=-0.21. annual income yields an increase in automobility of about 79 km in the short run, and 131 km in the long run (Table 7). As income is a key factor of household car ownership, income-elasticities of automobility should vary widely depending on whether the level of car equipment is held constant (as here) or not. In [START_REF] Hensher | The demand for vehicle use in the urban household sector: theory and empirical evidence[END_REF], the income-elasticities of car use for households living in the Sydney urban area ranged from +0.05 to +0.14 according to the level of car ownership of households. Cited for comparison in their study, [START_REF] Greene | The Influence of the Price of Gasoline on Vehicle Use in Multi-vehicle Households[END_REF] and [START_REF] Mannering | A Dynamic Empirical Analysis of Household Vehicle Ownership and Utilization[END_REF] respectively estimated this elasticity at +0.13 and +0.11 for the United States. Therefore, our results (Table 6) do not differ much from those got by these authors.

Conclusions

This article focuses on modelling the annual mileage covered by households with their personal cars, that is, their "automobility". It sheds new light on the car dependence issue. The results that are presented put an emphasis on the rational addiction model of [START_REF] Becker | An empirical analysis of cigarette addiction[END_REF], which had not been applied on automobile data so far. Original in our context, this model lives up to expectations, when applied to describe the empirical automobility behaviour of French households in 2000. Indeed, the assumption of addiction to car use is not statistically rejected. Therefore, the assertions made in [START_REF] Dupuy | From the "magic circle" to "automobile dependence": measurements and political implications[END_REF] and [START_REF] Wickham | Scenarios for a sustainable society: Car transport and the sociology of embedded technologies[END_REF] (mentioned in the introduction) are reinforced by the microeconomic point of view. Then, the rational version of the addiction model proves to be more relevant than the myopic version. Indeed, both past and future household automobility are significant to explain the current automobility in the rational addiction model. Moreover, the inter-temporal rate of substitution is valuated at 17 per cent, which is a plausible value. Such results show that the household behaviour is consistent with a theoretical inter-temporal optimization scheme. This conclusion stands our work apart from earlier dynamic studies, which may have missed an important point in explaining car use demand. Indeed, models based on a first-order autoregressive specification are useful to derive short and long-run elasticities, but they also require the individuals to be myopic as it regards the future. Nonetheless, the myopic model is also reported for comparison with other studies. In France, the kilometric-price and income elasticities of household automobility derived from the rational addiction model are in accordance with expectations. The respective estimations are -0.23 and +0.10 for the short run, -0.37 and +0.16 for the long run. These figures do not diverge from existing results that were reported in the literature. According to the type, petrol cars are more sensible than diesel cars to a change in fuel price. For diesel cars, the elasticity of annual mileage with respect to diesel-oil price is estimated at -0.13 at short run, while that for petrol cars with respect to premium-petrol price is measured at -0.32. The long-run elasticities are about 1.6 times higher. Our results are sensible whatever they are computed on the basis of the myopic addiction model or the rational addiction model. As it regards the latter approach, it strengthens our recommendation to use it when data are available. There are several ways to improve our approach. First, it would be interesting to demonstrate the model using disaggregated data over a longer period of time. Second, it is assumed in the rational addiction model that the future level of automobility is known by the agents, although there is some kind of uncertainty about future. It would therefore be relevant to test for different specifications that pertain to expectations about future. Also, it would be interesting to compare our study with any other approach that incorporates the fact that decision-makers could be forward-looking agents.

  ), or the inner suburbs of Paris ( 2 it R = ), or the outer suburbs of Paris ( 3 it R = ), or the Provinces ( 4 it R =

  estimating the model by means of ordinary least squares would produce biased estimates.

  our case, C drop to zero under

  evaluated at the average cars in 2000, using the estimates of the rational addiction model and the price for fuels in 2000 (see footnote 8). 95% confidence intervals in brackets.

  The elasticities in the rational addiction model can be expressed as functions of these characteristic roots. Let k E be the rd k continuous variable of E . The short and long-run elasticities of the demand for C with respect to k E , valuated at the sample averages k E and C , and respectively denoted /

	1 ϕ -measure the effect on it 2 C induced by a shock on	1 C + and it	1 C -it
	respectively. CT C Ek e	and / LT C Ek e	are given by

TABLE 1 :

 1 Descriptive statistics of household characteristics

	Year	1999	2000	2001
	Variable	Average Std Dev. Average Std Dev. Average Std Dev.
	# cars	1.18	0.82	1.19	0.80	1.20	0.82
	# adults	1.86	0.74	1.86	0.74	1.87	0.75
	# employed adults	1.00	0.85	0.98	0.85	0.98	0.86
	# women	0.98	0.48	0.98	0.48	0.98	0.48
	# children	0.52	0.93	0.50	0.91	0.49	0.90
	# licences	1.55	0.82	1.56	0.81	1.57	0.81
	Annual income (10 3 EUR)	23.08	13.68	23.82	13.95	27.70	14.37
	Age of the chief:
	<40	0.34	0.47	0.31	0.46	0.29	0.45
	]40 ; 65]	0.42	0.49	0.42	0.49	0.42	0.49
	>65	0.24	0.43	0.27	0.44	0.29	0.45
	Residential location:
	Paris-city	0.05	0.22	0.05	0.22	0.05	0.22
	Inner suburbs of Paris	0.06	0.24	0.06	0.24	0.06	0.24
	Outer suburbs of Paris	0.07	0.25	0.07	0.25	0.07	0.25
	The Provinces	0.82	0.38	0.82	0.38	0.82	0.38
	Automobility (km) (motorized HHs only)	19279	12759	18563	12313	18189	12576
	Automobility (km) (all HHs)	15610	13752	15193	13240	14826	13371
	Source: 1999-2001 French Car Fleet panel (3010 households).

TABLE 2 :

 2 Descriptive statistics of car characteristics the personal cars described by the households of the panel. Statistics for fuel-ambiguous cars (2% of the car sample for 1999, 1% for 2000 and 2001) not reported. Source: 1999-2001 French Car Fleet panel.

	Year	1999		2000		2001
		Average Std Dev. Average Std Dev. Average Std Dev.
	Repartition
	Diesel cars	0.34	0.23	0.36	0.23	0.38	0.24
	Petrol cars	0.64	0.23	0.63	0.23	0.61	0.24
	Mileage (km)
	Diesel cars	17786	9014	17085	9099	16187	8662
	Petrol cars	11001	6564	10412	6138	10157	6172
	Energy efficiency (L/(100 km))
	Diesel cars	6.82	1.29	6.80	1.49	6.74	1.35
	Petrol cars	7.74	1.63	7.72	1.58	7.67	1.56
	Age (in years)
	Diesel cars	5.19	4.01	5.52	4.28	5.48	4.38
	Petrol cars	7.22	5.70	7.50	5.86	7.56	5.99
	# observations	3552 cars	3576 cars	3605 cars
	Note: all

TABLE 3 :

 3 Evolution of the average fuel prices at filling stations in France (euros per litre)

	Year	1999	2000	2001
	Diesel-oil	0.6890	0.8461	0.7958
	Premium-petrol	0.9825	1.1380	1.0877
	Source: calculations from the yearbooks of the French professional comity of oils (Comité Professionel Des
		Pétroles).

TABLE 4 :

 4 Selection model estimates

	Variables	coefficient	t-stat
	Annual income (10 3 EUR)	1.42	4.91
	Age of the HH chief (reference: [40 ; 65[)
	[18 ; 40[	-0.10	-0.80
	≥ 65	0.22	2.19
	HH location (reference: Paris city)
	Inner suburbs of Paris	0.68	3.81
	Outer suburbs of Paris	1.64	8.10
	The Provinces	1.65	11.13
	# driving licence owners	1.50	20.22
	# adults (except the head)	0.07	0.91
	# employed persons	0.11	1.43
	# women	-0.36	-4.05
	# children	0.15	2.44
	Intercept	-2.63	-12.99
	Note: estimation for wave 2000 of the panel, 3010 households. Dependant variable:	1 Y = if the household i
	owned at least one car (2631 cases), 0 otherwise (379 cases). Probit estimation.

TABLE 5 :

 5 Estimates of the myopic and rational addiction models for household automobility

								Variable	Myopic addiction model	Rational addiction model
	Past automobility	(	1 KM -t	)	0.307	(3.41)	0.346	(4.20)
	Future automobility	(	1 KM + t	)	-	0.295	(3.98)
	(Cumulated) kilometric price (	∑	v v KMP	)	-362.68	(-3.04)	-223.29	(-1.98)
	HH's annual income	87.37	(3.98)	46.37	(2.16)
	# cars ( ) NC for HH living in:
		Paris-city ( =1) R	8813.36	(4.45)	5514.74	(2.94)
		the inner suburbs of Paris ( =2) R	8540.45	(4.46)	4888.76	(2.63)
		the outer suburbs of Paris ( =3) R	10810.79	(5.50)	6310.44	(3.08)
		the Provinces ( =4) R	11038.71	(5.56)	6189.06	(2.96)
	(Cumulated) age of cars:
		diesel (	∑	v v AD	)	-24.75	(-0.34)	-12.39	(-0.18)
		petrol (	∑	v v AP	)	-191.68	(-4.04)	-104.05	(-2.28)
		imprecise (	∑	v AI	v	)	-57.04	(-0.46)	-30.26	(-0.28)
		age not given (	∑	v AN	v	)	-2439.82	(-3.20)	-1112.06	(-1.53)
	# driving licence owners	2138.45	(3.13)	1324.65	(2.16)
	# adults (except the HH chief)	510.27	(1.01)	409.90	(0.96)
	# employed persons	815.36	(2.07)	286.32	(0.76)
	# women								-2041.94	(-3.34)	-1158.26	(-2.16)
	# children							-174.637	(-0.61)	-297.28	(-1.16)
	Age of HH chief (ref : [18-40[)
		[40-65[								-1617.14	(-2.64)	-1376.77	(-2.51)
		≥ 65								-2648.16	(-3.68)	-1608.01	(-2.34)
	Selection correction factor ( ) λ	2430.34	(1.71)	1555.84	(1.26)
	Intercept								-1210.78	(-0.90)	-954.01	(-0.85)
	ρ									∞	17.04%
	R	²								0.87	0.90
	Fisher								F(19, 2611): 82.6	F(20, 2610): 115.5
	(p-value)								(0.00)	(0.00)
	Breush-Pagan			2 X (40): 777.0	2 X (40): 878.8
	(p-value)								(0.00)	(0.00)
	Hansen								2 X (21): 30.46	2 X (20): 23.22
	(p-value)								(0.08)	(0.28)

TABLE 8 :

 8 Elasticities of annual mileage of cars with respect to fuel prices

	Horizon	Car type (fuel used)	Elasticity

Remember the assumption in the rational addiction model that the inter-temporal rate of substitution is equal to the interest rate of the economy. Arbitrarily, one can consider that a plausible rate should range between 0% and

20%. For instance,Becker (1996, p.103) reported a rate of 15%, arguing it is a "quite reasonable" value.2 As income, used in[START_REF] Becker | An empirical analysis of cigarette addiction[END_REF] or[START_REF] Baltagi | Rational addiction to alcohol: panel data analysis of liquor consumption[END_REF] for example.

If the variable λ proves statistically non-significant, the correction can be ignored. This point will be more discussed in section 5.2.

"Car Fleet" is the literal translation from French of the original name of the survey, "Parc Auto".

Sofres is the acronym for "Société Française d'Études par Sondages".

This phenomenon has even continued on the same pace during the 2000's: this proportion was about 50% in 2005 and 55% in 2008.

-0.22. Because of a larger inter-temporal dimension in the rational addiction model, the longrun price-elasticity is higher (-0.37) than in the myopic model (-0.31). In level, a permanent raise of EUR 1 of the cost to achieve 100 km causes a decrease in the annual mileage per car of 380 km in the short run, and 623 km in the long run (Table 7). -623 [-1194;-138] Note: evaluation from the rational addiction model (Table 5). 95% confidence intervals in brackets.

Estimating price-sensitivities of car use has been the topic of many works. To cite a few, such figures were estimated and discussed in [START_REF] Hensher | The demand for vehicle use in the urban household sector: theory and empirical evidence[END_REF], [START_REF] Oum | Concepts of price elasticities of transport demand and recent empirical estimates: an interpretative survey[END_REF], [START_REF] Eltony | Transport gasoline demand in Canada[END_REF], [START_REF] Rouwendal | An economic analysis of fuel use per kilometre by private cars[END_REF], [START_REF] Johansson | Measuring the long run fuel demand of cars: separate estimations of vehicle stock, mean fuel intensity, and mean annual driving distance[END_REF], [START_REF] Berri | Etalement urbain et croissance de la circulation automobile des ménages en France. Projection aux horizons 2010 et[END_REF]. [START_REF] Graham | The Demand for Automobile Fuel: A Survey of Elasticities[END_REF] collected many existing results in the literature dealing with car use sensitivities to fuel price, which can be expected to be close to our kilometric priceelasticities. In their study, [START_REF] Goodwin | A review of new demand elasticities with special reference to short and long run effects of price changes[END_REF] is a referenced paper: based on four elasticities drawn from empirical works in the 1980's, the author reported an average fuel price sensitivity of automobile traffic of -0.16 for the short run, and -0.33 for the long run. [START_REF] Goodwin | Elasticities of road traffic and fuel consumption with respect to price and income: a review[END_REF] updated this result with empirical works published in the 1990's and in the early 2000's. The authors reported an average sensitivity of mileage with respect to fuel price of -0.10 for the short-term and -0.30 for the long-term. In our present application, the most relevant model to be compared with the literature is probably the myopic addiction model, since short and long-run elasticities are usually derived from first-order autoregressive specifications. In this model, the long-run kilometric price elasticity of household automobility, estimated at -0.31, emerges as a very plausible value.


Valérie Moreau-Villéger

Hervé Delingette

Maxime Sermesant

Hiroshi Ashikaga

Elliot Mcveigh

Nicholas Ayache

Building Maps of Local Apparent Conductivity of the Epicardium with a 2-D Electrophysiological Model of the Heart

Keywords: heart modeling, electrophysiology, inverse problem, data assimilation, parameter estimation, reaction-diffusion system

In this article we address the problem of estimating the parameters of an electrophysiological model of the heart from a set of electrical recordings. The chosen model is the reaction-diffusion model on the transmembrane potential proposed by Aliev and Panfilov. For this model, we estimate a local apparent 2D conductivity from a measured depolarization time distribution. First, we perform an initial adjustment including the choice of initial conditions and of a set of global parameters. We then propose a local estimation by minimizing the quadratic error between the depolarization time computed by the model and the measures. As a first step we address the problem on the epicardial surface in the case of an isotropic version of the Aliev and Panfilov model.

The minimization is performed using Brent method without computing the derivative of the error.

The feasibility of the approach is demonstrated on synthetic electrophysiological measurements. A proof of concept is obtained on real electrophysiological measures of normal and infarcted canine hearts.

I. Introduction

A. Motivations

Cardiac arrhythmias are the cause of considerable morbidity. Tachyarrhytmias can originate from ectopic foci of electrical depolarization or from abnormal conduction pathways in the myocardium.

The treatment of choice for patients with tachyarrhythmias is radio-frequency ablation, where the abnormal electrical focus or pathway is ablated by applying radio-frequency energy. For patients with ventricular asynchrony, the treatment of choice is biventricular pacing through a pacing device.

An electrophysiological study is performed prior to these interventions: an electrical measurement catheter is inserted into the appropriate chamber of the heart and the electrical activity of the endocardial surface is measured. When inspecting electrophysiological data, cardiologists often base their analysis on the depolarization and repolarization maps of the epicardium or endocardium and more specifically on isochrones associated with both instants. From those maps, expert eyes can detect different electrophysiological pathologies ranging from the existence of low conduction zones caused by infarcted tissue, to the occurrence of fibrillation caused by scrolling waves.

The aim of the research effort presented in this paper is to provide cardiologists with a quantitative and objective map of one or several parameters related to conduction pathologies for a better diagnosis and a better planning of therapies (early detection of damaged cardiac tissue, or a precise location of ectopic foci or reentry points leading to an optimal planning of radiofrequency ablation). To achieve this task we do not resort to a pure signal processing approach, where for instance conduction could be estimated from the distance between two isochrone curves.

Instead, we propose to compute this additional information by solving an inverse problem: finding the parameters of a cardiac electrophysiology model that can best explain electrophysiological observations (depolarization time). In this article, we present a first step toward this objective: a method to estimate a local apparent conductivity that will be defined in Section II-C with a 2-D electrophysiological model of the heart from surfacic electrophysiological measures. This work is a first step towards a 3-D study.

B. Context

In electrophysiology, there are usually three different types of electrical potentials that can be considered: the extracellular potential, the intracellular potential and the transmembrane potential. The cardiac cells are separated from their environment by a membrane, creating two DRAFT distinct electrical spaces, respectively the intra and the extracellular domains which have different ionic concentrations and therefore different potentials, the intra and the extracellular potentials.

This difference creates a transmembrane potential (TMP). When a cardiac cell receives a pacing signal, the modification of the transmembrane potential creates an action potential leading to the contraction of cardiac fibers. For this reason, the measurement and computation of the transmembrane potential is crucial to assess the cardiac function.

The state of the art in measuring the electrical activity of the heart is mainly electro-anatomical mapping [START_REF] Rudy | Inverse computation of epicardial potentials[END_REF], catheter based measures like the EnSite system (Endocardial Solutions) or contact electrodes on epicardium or endocardium such as the studies presented in [START_REF] Mcculloch | Experimental and computational modeling of cardiac electromechanical coupling[END_REF], [START_REF] Faris | A novel technique for cardiac electromechanical mapping with MRI tagging and an epicardial electrode sock[END_REF], [START_REF] Rhode | XMR guided cardiac electrophysiology study and radio frequency ablation[END_REF]. In this paper, the measures were performed using an epicardial sock on several dog hearts. This procedure is very invasive but the electrodes are very accurately located and we show results on normal and infarcted hearts. The infarcted regions have different conduction patterns and therefore we expect the parameters of the model to be different in these regions.

A wide range of models of the electrical activity of the heart [START_REF] Noble | Models of cardiac ventricular action potentials: iterative interaction between experiment and simulation[END_REF], [START_REF] Noble | Modeling the heart from genes to cells to the whole organ[END_REF], [START_REF] Belik | Computational methods for cardiac electrophysiology[END_REF], [START_REF] Krinski | Cardiac muscle models[END_REF], [START_REF] Hunter | The IUPS physiome project: a framework for computational physiology[END_REF] has been developed from very accurate cellular models like Luo and Rudy models [START_REF] Luo | A model of the ventricular cardiac action potential. Depolarization, repolarization and their interaction[END_REF], [START_REF]A dynamic model of the cardiac ventricular action potential. simulations of ionic currents and concentration changes[END_REF] to phenomenologic efficient models [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF], [START_REF] Fenton | Vortex dynamics in three-dimensional continuous myocardium with fiber rotationfilament instability and fibrillation[END_REF], [START_REF] Bernus | A computationnally efficient electrophysiological model of human ventricle cells[END_REF], [START_REF] Tusscher | A model for human ventricular tissue[END_REF]. Accurate models such as bidomain models or Luo-Rudy models provide excellent insight into the physiological phenomena creating the electrical activity of the heart [START_REF] Colli-Franzone | Monodomain simulations of excitation and recovery in cardiac blocks with intramural heterogeneity[END_REF] but are probably too sophisticated for our inverse problem. Actually, these models are designed to capture very subtle modifications in the shape of the transmembrane potential [START_REF] Clayton | Propagation of normal beats and re-entry in a computational model of ventricular cardiac tissue with regional differences in action potential shape and duration[END_REF] whereas we only measure here a depolarization time distribution. Moreover subtle models contain too many parameters to estimate with respect to the number of available measures. The FitzHugh-Nagumo model [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] fits our expectations and allows to perform reasonably fast computations of the transmembrane potential propagation. Aliev and Panfilov have developed a modified version of the FitzHugh-Nagumo equations suited to the cardiac transmembrane potential [START_REF] Aliev | A simple two-variables model of cardiac excitation[END_REF].

In order to compare the transmembrane potential (TMP) computed by the model with the measured extracellular potential, we consider depolarization time. The depolarization time is the moment when a cardiac cell gets activated, leading to fiber contraction. The repolarization time is the moment when a cardiac cell returns to its rest potential, leading to fiber relaxation. The depolarization time can be computed from the extracellular potential measures, as explained in [START_REF] Sermesant | Preliminary validation using in vivo measures of a macroscopic electrical model of the heart[END_REF], by taking the instant when the time derivative is the most negative (Figure 1). From the transmembrane potential point of view, the depolarization time at a given point is the first instant when the transmembrane potential is above a threshold (Figure 1). As it is rather tedious to extract DRAFT repolarisation time from measured data, we based our study on the depolarization time only.

Estimating parameters from patient specific data can be addressed as a data assimilation problem. Data assimilation methods [START_REF] Van Der Merwe | The square root unscented Kalman filter for state and parameters estimation[END_REF], [START_REF] Evensen | Parameter estimation in dynamical models[END_REF] improve dynamical models by combining them with real observations. This combination is done via an appropriate minimization of the difference between observations and simulated data. Data assimilation techniques are very popular in meteorology and oceanography. They were recently introduced in biomechanical studies of the heart as in [START_REF] Shi | Stochastic finite element framework for cardiac kinematics function and material property analysis[END_REF] where an extended Kalman filter is used to estimate mechanical parameters of the heart from a real heart motion. None of the classical methods of minimization used in data assimilation are truly suited for the model and the measures of our problem, we thus propose a specific minimization method.

In Section II, we first describe the methods of this study. First we detail the model and the details of the simulation of the electrical wave based on this model. We then specify initial adjustment including the choice of initial conditions and of a set of global parameters before detailing the local apparent conductivity estimation. In Section III, we present the results of this procedure. We first validate the global adjustment and the local estimation on simulated measures. We then show the results of the estimation on a normal and on an infarcted heart. In Section IV, we discuss these results. Finally in Section V, we summarize and present the perspectives of this work.

II. Methods

The estimation method that we propose is based on the comparison between simulated and measured depolarization time on the epicardial surface. In this section we first present the Aliev-Panfilov model, followed by a procedure to get initial conditions and a set of parameters that produce a simulation that globally fits the measures. Finally we present a local parameter DRAFT estimation algorithm. This local estimation is performed using Brent method to minimize the quadratic error between measured and simulated depolarization time.

A. Model description

We chose the following formulation for the model of Aliev and Panfilov.

ε∂ t u = ε 2 div (D∇(u)) + ku(1 -u)(u -a) -uz (1.a) ∂ t z = -(ku(u -a -1) + z) (1.b) (1)
In this formulation, u is a normalized transmembrane potential (between 0 and 1) and z is a variable modeling the repolarization, k controls the repolarization, ε controls the coupling between the transmembrane potential and the repolarization variable z and a controls the reaction phenomenon.

A 3D anisotropic model based on the Aliev-Panfilov system was developed in the context of the ICEMA 1 collaborative research action [START_REF] Ayache | Towards model-based estimation of the cardiac electro-mechanical activity from ECG signals and ultrasound images[END_REF], [START_REF] Sermesant | Biomechanical model construction from different modalities: Application to cardiac images[END_REF], [START_REF] Chapelle | A physiologically-based model for the active cardiac muscle contraction[END_REF], [START_REF] Sermesant | Preliminary validation using in vivo measures of a macroscopic electrical model of the heart[END_REF].

The electrophysiological measures are usually available on the endocardium or the epicardium.

A first essential stage before addressing the 3D problem is to tackle a simplified and tractable problem by considering a surface model. Thus we perform a simulation based upon the Aliev and Panfilov model on a surface triangulation S with a set V of N vertices and a set T of L triangles.

The surface meshes used in this study are built from the locations of the electrodes of an epicardial sock.

We use a further simplification by considering an isotropic propagation i.e. D = d diag(1, 1, 1)

in system [START_REF] Rudy | Inverse computation of epicardial potentials[END_REF], where the diffusion coefficient d is proportional to a conductivity. The system (1) is scaled spatially to the maximum dimension of the triangulated mesh and temporally so that the action potential duration is around 0.3 s (τ = 0.26t, with τ the real time and t the normalized time). The temporal integration of the system (1) is done with an explicit Euler scheme. The spatial integration is performed with the finite elements method with linear triangular elements.

The numerical issues and the implementation are described in [START_REF] Moreau-Villéger | Global and local parameter estimation of a model of the electrical activity of the heart[END_REF]. Once a simulation of the entire cardiac cycle is performed, based on the Aliev-Panfilov model, the depolarization time distribution is computed using the following equation:

t s = min{t|u(t) = 0.5} (2)
1 ICEMA/ICEMA-2 : http://www-rocq.inria.fr/sosso/icema2/icema2.html As stated by [START_REF] Keener | A geometrical theory for spiral waves in excitable media[END_REF], the velocity of the depolarization wave on a 1D domain can be expressed as follows.

c = √ 2kd(0.5 -a) (3)
In 2D, the velocity of the depolarization wave is not constant in space. At each point in the mesh, it is equal to the velocity in 1D (Equation ( 3)) minus a term proportional to the curvature of the front [START_REF] Keener | A geometrical theory for spiral waves in excitable media[END_REF]. As we only need a global estimate of the propagation velocity on a surface, we neglect, as a first approximation, the front curvature and simply approximate the velocity c of the depolarization wave by its expression in Equation ( 3).

The depolarization velocity c can also be estimated from the gradient of the measured depolarization time on the surface, ∇ x t : 1/c = ∇ x t . The gradient ∇ x t on the surface can be computed as described in Appendix. Then, we compute a median value of the gradient ∇ x t over the whole mesh and using Equation (3) we get the following approximation: median (

∇ x t ) -1 ≈ √ 2kd(0.5 -a).
From this equation, if all parameters but k are fixed, it is possible to compute k directly:

k = 1 2d median ( ∇ x t ) -1 0.5 -a 2
But by doing so, we compare a theoretical 1D velocity and an apparent velocity computed on a 2D surface and this leads to a biased estimation. Therefore, we also compute a velocity estimated from a first simulation on the same mesh as the one used for the measures. As the velocity c is proportional to 1/ √ k, a ratio between measured c m and simulated c s propagation velocity can be DRAFT computed as follows.

median

∇ x t m median ∇ x t s = c s c m ≈ √ k s √ k m . (4)
The measured and simulated depolarization time distributions are denoted by t m and t s respectively. k s is the initial value for the parameter k in the initial simulation and k m is the new adjusted value. k m can be computed as follows.

k m = k s median ∇ x t s median ∇ x t m 2 (5)
With this new value for the parameter k, we get a complete set of global parameters that leads to a new simulation based on the model. The result of this simulation is the initialization of the local estimation that is described in the next section.

C. Local estimation of apparent conductivity

Once the simulated depolarization time map globally fits the measured one, a local adjustment of the model is possible. Depolarization time is an implicit function of the action potential whereas classical data assimilation methods generally require a linear function or at least an explicit function of the results of the model between the observations and the state variables. For that reason, we propose a specific estimation algorithm.

In the case of an infarct, all parameters a, k and d are modified in the infarcted region. The depolarization time distribution depends on the product √ kd(0.5 -a) and thus only allows to estimate one parameter. We choose the diffusion coefficient d as the spatially varying parameter.

Ischemic or infarcted regions exhibit a decrease in electrical conduction so that the velocity of the depolarization wave is smaller in this region and therefore should have a lower diffusion coefficient d.

Due to the fact that parameters a and k also vary in ischemic regions, the diffusion coefficient d also reflects the composite variations of parameters k and a. With the global adjustment of k proposed in section II-B, the diffusion coefficient d can no longer be related to the electrical conductivity.

Therefore, we call d the apparent conductivity (AC).

In the discretized model [START_REF] Moreau-Villéger | Global and local parameter estimation of a model of the electrical activity of the heart[END_REF], an AC value is assigned to each triangle. Consequently, we look for an AC map (d) = (d j ) 0≤j≤L-1 , where L is the number of triangles in the triangulation. This AC map should minimize

C(d) = v∈V (t m v -t v (d 0 , . . . , d L-1 )
) 2 where V is the set of the vertices in the triangulation, t m v is the measured depolarization time at vertex v and t v (d 0 , . . . , d L-1 ) the depolarization time at vertex v resulting from a simulation with the conductivities (d 0 , . . . , d L-1 ).

DRAFT

In order to have a robust estimate of the AC, we split the heart surface into different connected regions and estimate one AC value d R for each region R. Let (R k ) 0≤k≤K-1 be a partition of the surface in K regions. For each region R k , d j = d Rk for all j such that the j th triangle of the surface belongs to R k . Then, the new minimization problem is to find (d) = (d Rk ) 0≤k≤K-1 that minimizes

C(d) = v∈V (t m v -t v (d R0 , . . . , d RK-1 )) 2
In order to simplify the minimization of C, we want to use the causality of the electrical wave propagation to estimate the AC one region after another, following the order of depolarization.

To check this assumption, we studied the influence of the apparent conductivity d R of one region R on the depolarization time distribution by computing a rough estimation of ∂t/∂d R based on centered finite differences. For each vertex v ∈ V,

∂t v ∂d R ≈ t v (d R + ∆d) -t v (d R -∆d) 2∆d .
In order to compute t v (d R + ∆d) and t v (d R -∆d), we have to set the AC of region R to d R + ∆d and then to d R -∆d and perform one simulation for each of these values. We noticed that, as expected, the AC d R mainly influences regions having depolarization time equal to or greater than the depolarization time of region R.

As a consequence, we transform a K-dimensional minimization problem to K successive one dimensional minimization problems:

C(d R ) = v∈V (t m v -t v (d R )) 2 (6)
During the estimation of the AC value d R , the conductivities of other regions remain constant.

As there are enough vertices in a region to provide a robust estimate, we simplify the criterion C(d R ) by taking into account only the vertices of the region R: three points on the curve, there exists a unique parabola f (x) = αx 2 + βx + γ going through these points reaching its extremum at a point d x such that We construct an iterative process which is a simplified version of Brent's method [START_REF] Press | Minimization or Maximization of Functions[END_REF], to find value k m = 12.51. Figure 3 displays the error between the measured depolarization time and the simulated depolarization time before and after the estimation of k.

C(d R ) = v∈V,v∈R (t m v -t v (d R )) 2 (7)
d x = d b - 1 2 (d b -d a ) 2 (C(d b ) -C(d c )) -(d b -d c ) 2 (C(d b ) -C(d a )) (d b -d a )(C(d b ) -C(d c )) -(d b -d c )(C(d b ) -C(d a )) . (8)
We then test the local estimation algorithm on simplified simulated data. We still use the surface mesh of Figure 2. All the parameters but d are fixed: ε = 0.1, a = 0.15, k = 12 and dt = 5.10 -4 . We choose only one region R where the AC value is modified from the standard value 1.0 to a value d m R . On the rest of the mesh, the AC is equal to 1. We perform a new simulation with the modified AC and produce a map of depolarization time that we consider as input to the AC estimation algorithm. In this example, we suppose that we know the region R of modified conductivity.

We then present two examples, one with a low conductivity value and one with a large one. In this first example, the expected AC value is d m R = 0.5, we choose a very rough initial bracketing, d a = 0.1, d b = 2.0 and d c = 9.0. We perform the iterative process until the stopping criterion p is lower than 10 -3 . 33 iterations are necessary to obtain an estimation of d R = 0.501. Figure 4(a) shows the evolution of the estimated at each iteration of the minimization. The convergence is fast and stable. In a second example, the expected AC value is d m R = 5.0. We perform the iterative process until the stopping criterion is less than 0.01. By choosing the initial bracketing, d a = 2.0, d b = 6.5 and d c = 9.0, we assume that the AC d m R is greater than 2.0. With these parameters, 17 iterations are necessary to reach convergence and obtain d R = 5.02. As the function C(d R ) is very flat around its minimum, the convergence is slower with the parameters of the first example, but the estimated AC is the same. In this section, we explain in details the procedure for a normal heart. The in vivo measures used in this section were acquired on an adult male mongrel dog using a multi-electrode epicardial sock. The surgery, experimental layout and the data acquisition are described in [START_REF] Moreau-Villéger | Global and local parameter estimation of a model of the electrical activity of the heart[END_REF], [START_REF] Sermesant | Preliminary validation using in vivo measures of a macroscopic electrical model of the heart[END_REF]. This heart was artificially paced and the natural pacing was suppressed. For analysis of electrical activation, epicardial readings from each electrode were averaged over approximately 20 heartbeats.

The derivative of the voltage v(t) is computed with a five-point finite difference estimate. The depolarization time is chosen as the instant of the most negative derivative. The depolarization time distribution computed from the potentials recorded on the 128 electrodes of the sock is interpolated on a surface mesh consisting of 336 triangles and 192 vertices as displayed on Figure 5.

The pacing is artificial, consequently it is simple to select the initial pacing area as we show on Before the estimation, the mean error was 20.6 ms. After the automatic estimation, the mean error is 10 ms.

We now apply the proposed method to perform the local estimation of the AC. We first create a partition of the epicardium by splitting the epicardium into regions according to the isochrones of the depolarization time map. The precision of this partition is limited by the mesh resolution. In case this regions are too large we split them orthogonally to the isochrones. Finally, for this example we obtain 14 regions as shown in Figure 8.

We then iteratively estimate one AC value for each region. We sort the 14 regions of Figure 8 according to their depolarization time. The convergence on each region is presented in Figure 9 while Figure 10 shows different views of the obtained AC map.

The absolute error decreases significantly both after the global estimation and the local

C. Case of an infarcted heart

We also apply the AC estimation method on the case of an infarct on the anterior wall. This heart was artificially paced and the measures were obtained from an epicardial sock like in the normal case. The measured depolarization time distribution in this case is depicted in Figure 11(a).

The conduction appears to be slower in the infarcted region due to the local modification of the cells properties and of the conductivity. As a consequence, we expect the AC to be low in the

IV. Discussion

A. Performances and parameter choice.

The estimation of the AC for one region takes about 6 minutes on a dual processor Pentium The convergence of the local estimation algorithm is as stable for the real measures as for the simple synthetic example. The criterion C is approximated by a parabola, which exposes only one minimum, the global one. Therefore, we do not observe oscillations between several possible local minima. We did not experiment the drawbacks described in [START_REF] Press | Minimization or Maximization of Functions[END_REF]. Thus the simplified version of Brent minimization algorithm seems sufficient to minimize C(d).

This minimization method relies on a few parameters, the initial bracketing and the desired stopping criterion. A difference between two successive iterations lower than 0.01 is sufficient. In a first stage, we consider that no a priori information is available for the AC. A very rough initial bracketing is essential in order to assure that it contains the minimum of C. The AC should always be positive. Therefore, we take a small value for the lower bracket d a . Hence, the positivity of the AC d is easily imposed. There is no particular maximum limit for the AC. As all the regions belong to the same heart, it is realistic to take an AC always lower than 10.

An important stage of the estimation is the splitting of the epicardium in several regions where the estimated AC is constant. For the moment, this partition is tuned manually.

B. Expected accuracy of the apparent conductivity estimation

C. Discussion on the results

The global adjustment is satisfying and provides a good initialization to start the estimation of the local parameter. At this stage, we are more interested in the error on the depolarization time than on the value of k. The first reason is that k does not correspond to any measurable quantity that can be extracted from a set of measures. The second reason is that this global estimation is based on a model with spatially constant parameters whereas it may not be the case for simulated or real measures. The error is globally smaller after the estimation of k but the error is still large in regions far from the pacing region. Indeed, the estimation of the parameter k is not performed by controlling directly the depolarization time but by controlling the norm of the spatial gradient of the depolarization time derivative, that is the inverse of the velocity (Section II-B).

The estimated AC values are consistent with the measures. From Figure 9, we visually distinguish two types of regions. The regions of the first type have AC values around 0.9 and for a second class of regions, the AC is around 1.8. This reflects the asymmetry of the measures. The depolarization wave is faster on one side than on the other side of the epicardium. This may be due to the fact that we model the epicardium as homogeneous, without distinguishing the left and right ventricle nor taking into account the fiber directions. In agreement with the measures, our estimation algorithm provides us with asymmetric conductivities. On the upper left view of Figure 10, the green region corresponds to the pacing region. On the left, we estimate an AC close to 0.9, whereas on the right we estimate an AC up to 5.0. If we refer to Equation (3), it means that the velocity of the propagation is more than twice larger on one side than on the other ( 5/0.9 ≈ 2.34).

Let us consider more precisely the first region, which includes the pacing electrodes. In this region, we have a limitation with the global estimation because the depolarization appears to be faster than in other parts of the epicardium. The model with constant parameters is not able to capture this phenomenon. Fortunately with the local estimation, the AC value estimated in this region is greater than 1.0, producing a faster depolarization wave in the pacing region, in perfect agreement with the experiment.

From Figure 6, we notice a visual improvement of the depolarization time distribution. Compar-DRAFT ing these results with the measures (Figure 6(c)), we notice that the shape of the depolarization front resembles much the measures with the local adjustment. The quality of this estimation is also assessed by the visualization of the absolute error on the depolarization time in the epicardial surface. On Figure 7(b), the depolarization time error is small near the pacing region but becomes larger further away from this region. The final depolarization time error (Figure 7(c)) is not only smaller but also more homogeneous. Thus the agreement between the measures and the model is improved in the entire epicardial surface. In conclusion, the error has been considerably improved by the local estimation of the AC, even if we only added 14 degrees of freedom: 14 AC values for 14 regions.

D. Case with the infarct

On Figure 11(b), the infarcted region is displayed as bright circles. A large portion of the infarct is detected in the two regions with the lowest conductivity values, but we see that a part of the infarct is not detected as a low conductivity region. The heterogeneous infarct geometry in the heart wall can explain this observation: the infarct can be transmural (i.e. extending from the inner surface to the outer surface) or non-transmural (i.e. extending from the inner surface to somewhere in the wall), and when considering vertices in the mesh, where the infarct is non-transmural, electrical conductivity can be almost normal. In addition, a low conductivity is estimated in normal regions.

As seen in the first case, this may be due to the modeling of the epicardium as homogeneous medium. We are currently working on the inclusion of the fiber directions in this model.

V. Conclusions and perspectives

We addressed the problem of estimating a set of parameters for the transmembrane potential propagation modeled by Aliev and Panfilov from measured depolarization time. In order to evaluate the quality of our results, we used a criterion based on the difference of depolarization time between the model and the measures. We first presented a procedure to estimate globally a set of parameters so as to properly scale the electrical propagation. We used the theoretical properties of the Aliev and Panfilov equations and validated that approach on simulated depolarization time. We then presented a method to locally estimate the apparent conductivity (AC) region by region. This method was validated on simple simulated depolarization time. We successfully estimated global and local parameters of the model from in vivo measures of a canine heart. The simulation based on the model with these new values showed that the depolarization time error was significantly DRAFT decreased. Moreover, the AC values that we obtained are consistent with the measures even if a ground truth is not available. We also applied the estimation of the local AC to the case of an infarcted heart and found a strong correlation between parts of low AC values and the infarcted region.

In order to have a totally automatic process, we still have to find and build algorithms for the automatic partitioning of the epicardium. The proposed method could be easily extended to include more input data, like the isochrones of repolarization time, thus leading to the estimation of one additional spatially varying parameter (in this case the a parameter). In other words, we presented a general framework that is applicable to one or more macroscopic observations of cardiac electrophysiology. The model based approach that we propose could also be extended to a 3D model and could easily take into account fiber directions. As an isotropic model is too limited, we first plan to introduce anisotropy in the 2D model. The next step will be to estimate the parameters of a 3D heart model by establishing a correspondence between 2D measures and a 3D mesh. The proposed approach is not dependent on the model since it only uses simulations based on the forward model.

Thus, it can be adapted to more complex models that can reproduce specific pathologies. It can also be adapted to models computing depolarization time rather than transmembrane potential like eikonal models [START_REF] Colli-Franzone | Spreading of excitation in 3d models of anisotropic cardiac tissue :i. validation of the eikonal model[END_REF]. Finally a truly physiological validation would require to apply our method to a benchmark of pathological and normal measures analyzed by experts.

The electrophysiological measures used in this article are very rich but also very invasive.

Therefore, a long term challenge will be to estimate parameters of a full electromechanical model using both ECGs and displacement data, for example tagged MRI [START_REF] Axel | Tagged magnetic resonance imaging of the heart: a survey[END_REF]. This problem is more complex because the equations modeling the mechanical deformation are more unstable than for the electrophysiological models. First results concerning the estimation of mechanical parameters from MR images were obtained recently in [START_REF] Sermesant | Cardiac function estimation from MRI using a heart model and data assimilation: Advances and difficulties[END_REF].

Fig. 1 .

 1 Fig. 1. Comparing simulated transmembrane potential (left) and measured extracellular potential (right) on one electrode.

DRAFTB.

  scheme even for very coarse meshes associated to the epicardial sock. The parameters of the model a, k, or d can vary between individuals or species. With the knowledge of the depolarization time, we can only adjust one parameter. We choose in a first stage to estimate a global value for the parameter k from the depolarization time while standard values are assigned to the other parameters (a = 0.15 and d = 1.0). The diffusion coefficient d is locally estimated in a second stage (Section II-C).

  Many popular criterion minimization methods are based on the computation of the derivative of the criterion. In our case, it would involve the computation of the derivative of C and consequently of ∂t/∂d R . The computation of the derivative of the depolarization time at a given point with simple finite differences is very tedious and requires two simulation steps which are computationally expensive. More subtle methods like adjoint methods are not directly operative since the depolarization time is not an explicit function of the transmembrane potential. Therefore, we choose a minimization method that does not involve any derivative, an iterative inverse parabolic DRAFT interpolation like the Brent method[START_REF] Press | Minimization or Maximization of Functions[END_REF]. This very consistent method replaces the function to be minimized by a well-chosen parabola. The minimum of the function C is approximated by the easily and efficiently computed minimum of the parabola. Given (d a , C(d a )), (d b , C(d b )) and (d c , C(d c ))

  Then C(d x ) can be computed using a new simulation based upon the model of Aliev and Panfilov with this new value for d and compared with C(d a ), C(d b ) and C(d c ).

Fig. 2 .Fig. 3 .

 23 scheme converge. The depolarization time distribution of Figure 2 is now a set of simulated data. Applying the procedure described in Section II-B to these simulated data, we obtain a global value of k m = 12.51 starting at a crude initialization of k s = 8. The values of a and ε for the crude initialization were chosen equal to the parameters used to simulate the data and the AC d was chosen spatially constant and equal to 1.0. A simulation is then performed with this new k value and the depolarization time distribution is compared with the simulated measures. The depolarization time error drops from 25.3 ms with k s = 8 to 14.7 ms with the new estimated

DRAFT

Fig. 4 .

 4 Fig. 4. Evolution of the estimated apparent conductivity during the estimation. (a): d m R = 0.5, p = 10 -3 . (b): d m R = 5.0, p = 10 -2 .

Figure 4 (

 4 b) shows the iterative estimation of AC values. B. Results on in vivo measures. Normal case.

Figure 5 .

 5 Figure 5. All the parameters but k are fixed as explain in section II-B, a = 0.15, ε = 0.1, d = 1 and the time step is dt = 5.10 -4 . Applying the method presented in Section II-B to the data, we obtain a global value of k m = 25.2 starting from a crude initialization k s = 8.

Figure 6 (Fig. 5 .Fig. 6 .

 656 Figure 6(a) presents the depolarization time obtained after the global adjustment. The depolarization time error obtained before and after the automatic global estimation procedure are

Fig. 7 .Fig. 8 .

 78 Fig. 7. Absolute depolarization time error between measures and simulations before (a) and after (b) the global automatic estimation and finally after the local estimation (c)

Fig. 9 .

 9 Fig. 9. Convergence of the apparent conductivity for the 14 selected regions.

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Apparent conductivity map. Different views. The color represents the apparent conductivity value. From blue (small values) to red (high values).

Figure 7

 7 Figure 7 displays the absolute error on the depolarization time before and after the local estimation of the AC.

  infarcted region. The AC values are displayed on Figure 11(b), the bright circles correspond to the approximate location of the infarcted region. The depolarization time distribution computed from a simulation taking into account these values is displayed on Figure 11(c). In the infarcted region, the shape of this depolarization front reproduces the shape of the measured depolarization front (Figure 11(a)).

III 1 Fig. 12 .

 112 Fig. 12. (a) Depolarization time as a function of the apparent conductivity dR at the seven vertices of a region R. (b) the derivative ∂t/∂d of these curves

Figure 12 (

 12 Figure 12(a) presents an example of the depolarization time as a function of the AC d R at the seven vertices of a region R and Figure 12(b) presents the derivative ∂t/∂d of these curves. When the absolute value of the derivative of t according to d is small, it means that the depolarization time distribution does not change much upon a moderate variation of the AC. The larger is the derivative, the smaller the range of AC values that leads to approximatively the same depolarization

Fig. 13 .

 13 Fig. 13. Notations used on a triangle (P0P1P2).

Acknowledgments

This work was partially funded by the scientific direction of INRIA through the Cooperative Research Action ICEMA2 (http://www-rocq.inria.fr/sosso/icema2/icema2.html). This action involved several partners including the INRIA research groups CAIMAN, EPIDAURE, MACS and SOSSO and was coordinated by Frédérique Clément (SOSSO).

DRAFT

Appendix

On each triangle trngl = (P 0 P 1 P 2 ) of the triangulation, the gradient of t is

where A is the area of the triangle (P 0 P 1 P 2 ) and S i = l i n i , l i is the length of the segment P i+1 P i+2

and n i is the external normal to this segment (Figure 13). We denote SH(v) the shell of the vertex v, that is the set of triangles to which v belongs. Then for each vertex v ∈ V, the gradient ∇ x t is given by

where A(trngl) is the area of triangle trngl.