Louise Boisen

Ryun Kim

Anna Sofie Drasbek

Peter Pedersen

Kristensen

Peter Kristensen

Evaluation of endothelial cell culture as a model system of vascular ageing

Keywords: Ageing, Endothelial cells, senescence, tube formation, SA beta-galactosidase activity, lysosome staining, gamma-H2A.X staining

The purpose of this study was to evaluate the relevance of long-term endothelial cell culture as a model system of vascular ageing. Micro-and macrovascular endothelial cells were serially passaged until replicative senescence and their ability to form tube-like structures when cultured on Matrigel was assessed throughout their lifespan. For both cell types low passage cultures adopted a homogeneous cobblestone morphology, while senescent cultures were extremely heterogeneous.

Furthermore, both cell types showed a reduction in tube formation ability with in vitro ageing, which is in accordance with the reduction in angiogenic potential observed with ageing in vivo.

Examination of senescence associated β-galactosidase activity revealed an increased activity in cells forming tubes as compared to cells cultured on plastic, which could be attributed to an increased lysosomal content of cells undergoing tube formation. As this increased senescence associated βgalactosidase activity was unrelated to the replicative age of the cells, senescence associated βgalactosidase activity may not be a relevant senescence marker for differentiating endothelial cells.

The age-related reduction in tube formation ability suggested that long-term culture of endothelial cells may be a valid model system of vascular ageing, which makes it an ideal platform for high throughput screening of compounds influencing angiogenesis.

A C C E P T E D M A N U S C R I P T

Introduction

Common to all vessels of the vasculature is the endothelium, which is a single layered sheet of squamous, polarised cells that are primarily responsible for all signalling and transportation from blood to tissue and vice versa [START_REF] Cines | Endothelial cells in physiology and in the pathophysiology of vascular disorders[END_REF]. As the body ages, so does the vasculature and ageing of the vasculature is accompanied by a number of changes such as reduced vasoregulation (Taddei et al, 2001), delayed and/or altered angiogenesis [START_REF] Agah | Thrombospondin 2 levels are increased in aged mice: Consequences for cutaneous wound healing and angiogenesis[END_REF], Sadoun & Reed, 2003), reduced repair capacity and altered blood vessel composition [START_REF] Izzo | Aging and arterial structure-function relations[END_REF].

Different types of endothelial cells have been isolated and cultured to understand the biology and pathobiology of the vasculature and the angiogenic response. The most widely used type isolated from human umbilical vein (named HUVEC) show a characteristic cobblestone shape when cultured, although cell borders are indistinct and hard to visualise (Jaffe et al, 1973). Senescent HUVEC cultures have been described as very heterogeneous with greatly enlarged cells containing multiple cytoplasmic protrusions coupled with positive staining for senescence associated βgalactosidase (SA β-gal) activity (Kalashnik et al, 2000, Unterluggauer et al, 2003).

In vivo, angiogenesis is mediated by the microvasculature and the angiogenic response of human microvascular endothelial cells (HMVEC) has been suggested to be stronger than that of HUVEC because of differing expression of matrix metalloproteases (Jackson & Nguyen, 1997), However, HMVEC are generally more fastidious in their culture requirements than HUVEC making longterm culture more challenging.

Functional studies of endothelial cells have been carried out in culture, as endothelial cells cultured in the presence of an extracellular matrix (ECM) organise themselves into tube-like structures.

Depending on the conditions, the tube structures will form in two or three dimensions (Staton et al, 2009). This response is seen as a mimic of natural endothelial behaviour during angiogenesis, and as such has been used as a simple in vitro angiogenesis assay. As expected for endothelial cells, low

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

passage HUVEC form tube-like structures when cultured on top of gels formed by ECM proteins such as collagen-1 or Matrigel (Unger et al, 2002). However, such studies have not been reported

for micovascular endothelial cells, although these cells are directly involved in angiogenesis. While endothelial cells in vivo normally are quiescent, cultured endothelial cells are kept in an artificially activated state characterised by extensive proliferation. Therefore, cultured endothelial cells are not ideally suitable for the study of the quiescent state. In contrast, these cell cultures are in principle well suited for the study of activated endothelium e.g. the study of angiogenesis. Furthermore, as primary cells in culture eventually reach replicative senescence, functional studies of cultured endothelial cells could turn out to be an excellent model system of angiogenesis during ageing.

In this study, we have cultured both micro-and macrovascular endothelial cells until the onset of replicative senescence. Minor differences in morphology between the two cell types at corresponding ages were found, while the senescent cells were distinguishable from young cells on the basis of morphology and SA β-gal activity. Finally, studying the function of the activated endothelial cells, we found that in vitro tube formation of the cells diminished as the culture aged and, thus, the number of senescent cells increased. This observation suggests that long-term culture of endothelial cells may serve as a suitable model of vascular ageing. Intriguingly, a pronounced and increasing SA β-gal activity was observed in cells regardless of age when cultured on Matrigel.

However, this increase was not accompanied by a corresponding increase in γ-H2A.X positive cells after cultivation on Matrigel. This suggests that the increased SA β-gal activity did not reflect a similar increase in senescent cells and questions the validity of this particular marker in the evaluation of senescence in differentiated endothelial cells.

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Materials and Methods

Cell culture

HMVEC isolated from adult dermis as well as HUVEC were obtained from Clonetics Cell Systems (Lonza, Basel, Switzerland For each passaging the resulting PD estimation was added to the sum of PD from the previous passages to achieve the cumulative population doubling level (CPDL), which was plotted against time to obtain a growth curve.

Tube formation assay

Wells of 96-well (Corning, Lowell, MA) or 4-well (Nunc, Roskilde, Denmark) cell culture plates were coated with BD Matrix Growth Factor Reduced Matrigel™ (BD Biosciences, San Jose, CA) at 4°C. The gel was allowed to polymerise in the 37°C CO 2 -incubator for at least 30 minutes before

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

seeding of 20 000 (96-well) or 80 000 (4-well) cells on top of the gel. The cells were allowed to attach and differentiate overnight. Tubes were visually examined with phase contrast using either 2.5x or 10x objective and pictures were recorded digitally with a Zeiss Axiocam MRc camera (Carl Zeiss Micro Imaging, Göttingen, Germany) attached to the microscope. Due to Matrigel coating and light refraction through the gel and media, dark spots in the centres and sometimes edges of wells were visible in the pictures.

Senescence evaluation -SA β-gal activity assay

Cells were examined for senescence associated β-galactosidase activity, which is a widely accepted marker for cellular senescence [START_REF] Dimri | A biomarker that identifies senescent human cells in culture and in aging skin in vivo[END_REF]. Briefly, cells were washed twice with phosphate buffered saline (PBS) and fixed with 2% formaldehyde, 0.2% glutaraldehyde solution in PBS for 3-5 min at room temperature (RT). The fixed cells were washed twice with PBS before addition of 1 mg/ml X-gal in NN-dimethylformamide, 2 mM MgCl 2 , 150 mM NaCl, 5 mM potassium ferricyanide, 5 mM potassium ferrocyanide, 40 mM citric acid/Na phosphate buffer, pH 6.0. Cells were incubated with the staining solution overnight at 37°C and washed twice with H 2 O and allowed to dry before visual examination in an inverted light microscope. Cells forming tubes on Matrigel were stained using the same procedure, although the stained tubes were not dried prior to examination.

Lysosome detection

The lysosomal content of the cells was evaluated using a fluorescent LysoTracker probe (Invitrogen, Carlsbad, CA). HUVEC NS2 cells were seeded on plastic and Matrigel in 24-well plates (Corning) at a density of 40 000 cells per well and allowed to attach/form tubes overnight.

The following day the growth medium was replaced with fresh, warm growth medium containing 100 nM LysoTracker Red DND-99, and the cells were incubated for 2 hours in the CO 2 -incubator.

The LysoTracker solution was replaced with normal growth medium and the cells examined for red

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

fluorescence under an inverted fluorescence microscope (Leica DMI 3000B (Wetzlar, Germany)

equipped with an Olympus DP72 CCD camera (Olympus, Ballerup, Denmark).

Kinetic study of LysoTracker and SA β-gal staining

To study the staining of lysosomes with LysoTracker Red DND-99 and SA β-gal activity, low and high passage HUVEC NS2 cells (40 000 cells/well) were plated in 24-well plates with or without prior coating with Matrigel and assayed 4, 8, 12, and 16 hours after plating. At each time point, live cells were stained with LysoTracker for 30 min before fluorescent pictures were obtained with the same exposure conditions (200 ms exposure) at all time points using an inverted fluorescence microscope. Simultaneously, the wells for SA β-gal activity visualization were fixed at each time point and incubated with X-gal solution for 16 hours before acquisition of pictures using differential interference contrast (DIC).

Immunocytochemical staining

Low and high passage HUVEC NS2 cells were seeded in 24-well plates with or without prior coating with Matrigel and allowed to grow for 16 hours before fixation in 2% paraformaldehyde diluted in PBS for 15 min. The cells were permeabilized with 0.1% Trion-X-100 for 5 min prior to blocking in 5% FCS for 1 hour. After washing, the wells were incubated with primary antibody, rabbit-anti-γ-H2A.X (Abcam, Cambridge, UK) for 1 hour followed by incubation for 45 min with goat-anti-rabbit-Alexa546 (Invitrogen). Stained cells were overlaid with VectaShield antifade with DAPI (Vector Laboratories, Burlingame, CA) before visualization using an inverted fluorescence microscope equipped with filter sets for DAPI and red fluorescence.

Cellular expression of endothelial cell markers

Cells were analysed for expression of endothelial cell markers von Willebrand Factor (vWF) and

Results

Cell cultures

The human microvascular cell line, HMVEC, was continuously passaged until senescence, which ensued after approximately 40 PD (Figure 1A). The growth curve showed the expected pattern of rapid growth for the first 35 PD, or about 40 days, where after growth slowed down over a period of about 20 days and eventually subsided, creating a plateau. As the growth curve also shows, the nearly non-growing culture was rather stable with no significant loss of cells over time. Low passage HMVEC cultures were very homogeneous with small, spindle-shaped cells growing in monolayers. When cells approached confluence they adopted the typical cobblestone morphology of endothelial cells (not shown). HMVEC appeared reasonably homogeneous with no significant SA β-gal activity for more than half of their lifespan (Figure 1E). The cellular rate of proliferation started slowing down around 30 CPDL (Figure 1A), which coincided with the appearance of enlarged and oddly shaped cells. The proportion of senescent cells increased in the last quarter of the culture's lifespan (Figure 1F) until the very end, where the culture was dominated by a heterogeneous SA β-gal positive population of greatly enlarged cells with multiple cytoplasmic protrusions (Figure 1D andG). The cells stained positively for endothelial markers, vWF and CD31, throughout their lifespan (not shown).

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

The long-term culture of HUVEC largely followed the pattern seen for HMVEC culture with the cells undergoing replicative senescence at about 36-38 CPDL after approximately 80 days in culture (not shown). Low passage HUVEC exhibited the characteristic cobblestone morphology described for endothelial cells, with the majority of cells being rather small and slightly spindle shaped. As the cultures aged, they gradually became more heterogeneous with oddly shaped cells starting to appear. Generally, cells of abnormal shape and size stained positive for SA β-gal activity, which resulted in a gradual increase in SA β-gal positive cells during ageing (not shown) and as seen in HMVEC cultures, the HUVEC culture eventually became dominated by greatly enlarged SA β-gal positive cells with multiple cytoplasmic protrusions when growth subsided (not shown). In contrast to the HMVEC cultures, there was a substantial cell loss in the last few passages, where the cells detached from the flask and died.

Tube formation assay

The tube formation assay was used to test the function of the activated endothelial cells throughout their lifespan. At several time points HMVEC were assessed for their tube formation ability as shown in Figure 2. Low passage cells formed many short tubes and the overall picture resembles that of a honeycomb (Figure 2A). As the cells aged and entered the second half of their lifespan, the tubes became longer and the branch points fewer (Figure 2B-E). Finally, as the growth rate slowed down and the culture reached senescence, the cells eventually stopped forming tubes altogether (Figure 2E andF) while the cells appeared to form clusters, which may represent the node points of the tubes.

As for HMVEC, the tube formation ability of HUVEC was affected by cellular ageing. While low passage (not shown) and mid-life (Figure 3A) cells did form tubes when cultured on Matrigel, this ability appeared to diminish as the culture approached senescence, with senescent cells hardly forming tubes at all (Figure 3C). As seen in Figure 3, the ageing cells appeared to concentrate in

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

branch points rather than spread out and form tubes. This behaviour could reflect a reduced ability to migrate and/or differentiate with advancing age.

Correlation between tube formation ability and appearance of senescent cells in culture

To see if the appearance of senescent cells in culture could be definitely linked with reduced tube formation, the SA β-gal assay was attempted for HUVEC NS2 cells cultured on Matrigel. Cells cultured on Matrigel did stain positively for β-gal activity in an age-related manner, but, interestingly, stained cells appeared earlier and were much more prominent when cultured on Matrigel than when cultured on plastic (Figure 4). Because SA β-gal activity has been shown to be lysosomal β-gal activity, it seemed relevant to evaluate whether the increased SA β-gal activity observed in the tubes could be caused by an increased lysosomal content. To achieve this goal, cells seeded on Matrigel and cells seeded on plastic were incubated with a LysoTracker probe that accumulates and fluoresces in the acidic organelles of cells. As seen in Figure 5, the tubes were strongly fluorescent (Figure 5C) while fluorescence in the cells seeded on plastic generally was much weaker with only the occasional cell showing strong fluorescence (Figure 5D). The picture for cells assayed for SA β-gal activity was similar, with tubes being strongly positive (Figure 5E)

and cells seeded on plastic less positive (Figure 5F).

To investigate whether the increased LysoTracker/SA β-gal activity in tube forming cells was a result of differentiating endothelial cells, the LysoTracker fluorescence and SA β-gal activity were evaluated at different time points after seeding in both low passage and high passage HUVEC. As LysoTracker intensity was also evident in the cells seeded on plastic (Figure 6, I-L). In parallel, cells were assayed for SA β-gal activity where a similar but less pronounced increase in activity with time was observed in cells seeded on Matrigel (not shown).

Evaluation of endothelial senescence with DNA damage marker during differentiation

As the increased SA β-gal activity in endothelial cells cultured on Matrigel indicated that SA β-gal activity might not be a good marker for endothelial cell senescence, low and high passage HUVEC were analysed for the presence of a DNA damage response signified by the expression of γ-H2AX.

As seen in Figure 7, high passage cells showed punctuate nuclear staining with a γ-H2AX antibody at a similar level whether cultured on plastic or Matrigel (Figure 7D and H), while nuclear staining of low passage cells was extremely limited (Figure 7B andF).

Discussion

Endothelial cell culture has been used extensively to describe many of the changes happening with ageing of the vasculature. However, the application of a continuously dividing cell culture system, which eventually undergoes senescence, to the study of a cell type that in vivo is mostly quiescent has some obvious disadvantages. This study was undertaken to investigate the validity of the Hayflick model on a particular aspect of vascular ageing, namely the well-documented reduction in blood vessel formation happening with ageing.

To cover different parts of the vascular endothelium both macro-(HUVEC) and microvascular (HMVEC) primary cell cultures were used. As described, HUVEC is the most extensively studied endothelial cell type, primarily because they are relatively easy to isolate and culture. Therefore, this cell type was chosen rather than more recently available and not thoroughly characterised somatic macrovascular endothelial cell types, although they might have been more interesting from

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

a physiological point of view. The dermal HMVEC isolated from adults were chosen because they represent a cell type that is often activated and involved in angiogenesis.

Morphologically, HUVEC and HMVEC cultures were quite similar in young and middle age with slight variations in shape from the traditional cobblestone appearance to more of a spindle shape, although the degree of confluence in the cultures did influence the cell shape. As the cells aged, they started showing typical signs of senescence with changes in shapes and increased size as well as SA β-gal activity. In fact, the most prominent difference between senescent HUVEC and HMVEC was the observed quite substantial cell loss during the last few passages of HUVEC cultures. This might be due to apoptosis as senescent HUVEC has been reported to undergo extensive apoptosis [START_REF] Hampel | Differential regulation of apoptotic cell death in senescent human cells[END_REF]. In contrast, apoptosis did not seem to be prominent in the senescent HMVEC cultures, which could be maintained for several weeks with no obvious changes in cell density. This is in accordance with earlier reports of stable senescent HMVEC cultures (Vasile et al, 2001, Watanabe et al, 1997). These observations indicate that not all types of endothelial cells behave alike in culture and that some cell-type specific characteristics may be retained in culture making the use of endothelial cells isolated from different parts of the vasculature appropriate.

During culturing, the functionality of the endothelial cells was tested using the tube formation assay. Interestingly, cellular senescence appeared to have a profound effect on cellular behaviour as tube morphology was influenced by the age of the assayed cells. For both HMVEC and HUVEC, cellular ageing was accompanied by a gradual increase in the length of the tubes coupled with a decrease in the number of branch points. As the cells entered the senescence plateau, tube formation was nearly completely abolished with the cells apparently forming clusters. Thus, these results

suggest that in vitro ageing of endothelial cells causes a loss of differentiation potential as evidenced by the reduction in tube formation from senescent cultures. Interestingly, most reports on

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
long-term cultures of endothelial cells do not assess tube formation in any form, but the few reports available have concluded that senescence of HUVEC does indeed lead to loss of tube formation ability [START_REF] Ha | Outgrowing endothelial progenitor-derived cells display high sensitivity to angiogenesis modulators and delayed senescence[END_REF][START_REF] Yang | Human endothelial cell life extension by telomerase expression[END_REF]. Reports on long-term culture of HMVEC are scarce and apparently none of them have addressed tube formation ability, although migration of HMVEC was reportedly reduced with ageing (Reed et al, 2000).

In an attempt to further investigate the correlation between endothelial cell senescence and tube formation ability, HUVEC NS2 cells were assessed for their ability to form tubes throughout their lifespan. To visualize the proportion of senescent cells, the tubes were stained for SA β-gal activity, as were control cells seeded on plastic. The staining procedure revealed that endothelial cells seeded on Matrigel did give rise to age-related SA β-gal staining although the apparent proportion of stained cells on Matrigel was higher than the proportion of cells grown on plastic. To further investigate the consistent difference between cells grown on plastic and Matrigel, another staining molecule was used as the SA β-gal assay does not measure cellular age per se but rather the abundance of lysosomes in the cytoplasm of the cells (Kurz et al, 2000). Thus, the early appearance of SA β-gal staining could be due to an upregulation of the number of lysosomes during cell differentiation when cultured on Matrigel. By tracking the acidic organelles including the lysosomes with fluorescent LysoTracker probes, it was found that tube forming cells exhibited a much stronger level of fluorescence than individual cells seeded on plastic, which correlated well with the difference in SA β-gal staining intensity. Moreover, it was found that the fluorescence intensity as well as the SA β-gal activity increased with time after seeding on Matrigel but not on plastic, which supported the idea that endothelial cells upregulate the number of lysosomes during differentiation. The observed evolution in fluorescence over time was seen in both high passage and low passage cells (Figure 6A-H) supporting this theory of increasing numbers of lysosomes.

However, no articles have been published supporting the notion that endothelial cells cultured on

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

Matrigel should contain more lysosomes than endothelial cells in standard culture. Alternatively, the difference in staining may be explained by the actual differentiation process endothelial cells undergo during tube formation: endothelial cells undergoing tubulogenesis in vivo and in 3D culture systems form a lumen through the production and fusion of multiple acidic cytoplasmic vesicles (Lubarsky & Krasnow, 2003, Ulmasov et al, 2009). Whether endothelial tube formation in two dimensions is accompanied by lumen formation is questionable, but Matrigel does lead endothelial cells down the differentiation pathway (Staton et al, 2009), and as such, it is not inconceivable that the endothelial cells will form intracellular acidic vacuoles when cultured on Matrigel. The vacuoles that contribute to lumen formation are not lysosomes (Ulmasov et al, 2009) and presumably do not contain β-galactosidase. However, the fact that they are acidic might explain why endothelial cells cultured on Matrigel exhibited a higher degree of SA β-gal staining than expected for their age. If SA β-gal is in fact normal lysosomal β-gal, the enzymatic activity is maximal around pH 4. The staining procedure is performed at pH 6.0. However, if the endothelial cells contain a considerable amount of acidic vesicles, that might be enough to shift the pH of the buffer sufficiently low as to allow a greater enzymatic activity and consequently a more intense staining. In accordance with this, the LysoTracker probe does not fluoresce in lysosomes specifically, but rather in acidic vacuoles in general.

The observation that SA β-gal activity increased in cells cultured on Matrigel regardless of age naturally questions the validity of this presumed marker of cellular senescence in assessing the cellular age of differentiating endothelial cells. In order to confirm that the observed increase was unrelated to the replicative age of the cells, cells cultured on plastic and on Matrigel were analysed for the expression of γ-H2A.X, which is a component of the DNA damage response and has been shown to be expressed in the nucleus of senescent cells (Wang et al, 2009). As expected, γ-H2A.X staining was found in the nuclei of high passage cells with a similar frequency regardless of the

A C C E P T E D M A N U S C R I P T

ACCEPTED MANUSCRIPT

culture conditions (Figure 7). In contrast, in low passage cells the staining was negligible both on plastic and Matrigel as would be expected based on the replicative age of the seeded cells.

Consequently, a DNA damage response does not appear to be initiated during differentiation on Matrigel. Previously, various studies have questioned the relevance of SA β-gal activity as a relevant marker of cellular senescence. [START_REF] Yang | The limitations and validities of senescence associated-β-galactosidase activity as an aging marker for human foreskin fibroblast Hs68 cells[END_REF] found that SA β-gal activity in human foreskin fibroblasts increased not only with cellular age, but also in confluent quiescent cells, in serum-starved cells and in cells treated with H 2 O 2 . Additionally, a number of studies have failed to find a correlation between the intensity of SA β-gal staining and donor age in vivo [START_REF] Cristofalo | SA β Gal staining: Biomarker or delusion[END_REF]. In conjunction with these observations and based on the presented data, SA β-gal activity as a measure of cellular senescence is questionable when assaying differentiating endothelial cells.

In conclusion, from the presented results it appears possible that senescence of the endothelium contributes to the age-related reduction in angiogenesis, and that long-term endothelial cell culture can serve as an appropriate, albeit simple, model of the reduction of angiogenic potential occurring during vascular ageing. This model will enable the high throughput screening of compounds and antibodies having an effect on angiogenesis during ageing. Furthermore, these results indicate that endothelial differentiation when cultured on Matrigel is accompanied by an increase in the number of lysosomes as illustrated by the LysoTracker probe and SA β-gal activity.

Additionally, this observation highlights that SA β-gal activity may not be particularly relevant as a marker of senescence during endothelial differentiation and should in any case not be used as the sole marker when evaluating endothelial senescence but should be confirmed by independent

  platelet endothelial cell adhesion molecule/CD31. Cells were washed trice in cold PBS before fixation for 2 min at -20°C with a methanol/acetone mixture. After three additional washes, the A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT culture flasks were blocked for 1 hr RT with 1% BSA. Cells were incubated with primary antibodies, rabbit anti-vWF (Abcam, Cambridge, UK) and mouse anti-CD31 (Abcam) shaking overnight at 4 °C. The cells were incubated for 2½-3 hrs with fluorescein isothiocyanate (FITC)conjugated anti-rabbit antibody (DAKO, Glostrup, Denmark) and Rhodopsin-conjugated antimouse antibody (Invitrogen) shaking in the dark before counterstaining with DAPI (Invitrogen) and examination for green, red and blue fluorescence under an inverted fluorescence microscope.

Figure 6

 6 Figure6shows, there was a clear increase in LysoTracker fluorescence with advancing time in the

Figure 1 .Figure 2 .

 12 Figure 1. Cultivation of primary human microvascular endothelial cells (HMVEC). A: growth curve showing the cumulative population doubling level (CPDL) during culturing from young to senescent. B-D: cellular morphology changes with advancing cellular age from small cobblestone shaped young cells to large cells with multiple protrusions at old age. E-G: the activity of SA β-gal increases with advancing cellular age. B: CPDL 22, C: CPDL 34, D: CPDL 39, E: CPDL 22, F: CDPL 35, G: CPDL 39. Scale bar 100 µm.

  ). HUVEC NS2 were kindly provided by Professor Pidder Jansen-Dürr of the Austrian Academy of Sciences in Innsbruck (Austria).

	HMVEC were cultured in Endothelial Cell Basal Medium-2 (Lonza) supplemented with EGM-2
	singlequots (Lonza) containing h-EGF, VEGF, h-FGF-B, R3-IGF-1, hydrocortisone, and FBS (2%
	final concentration). HUVEC were cultured in Endothelial Cell Basal Medium (Lonza)
	supplemented with EGM Singlequots (Lonza) containing h-EGF, Hydrocortisone, GA-1000, bovine
	brain extract (BBE), FBS (2% final concentration). HUVEC NS2 were cultured in EGM in cell
	culture flasks coated with 0.2% gelatine in EGM.
	The cells were grown in an atmosphere of 5% CO 2 at 37°C and subcultured by trypsinisation with
	trypsin-EDTA (Lonza). Cells were passaged such that the monolayers never exceeded 75-90%
	confluency. The cells were propagated until senescence and cell numbers were determined when
	subcultured by passing a 1/50 dilution of cells through a Coulter Counter. Population doublings
	(PD) were estimated using the following equation:
	PD	=	log	10	( #	cells	harvested	at	end	of	passage	) log log 2 10 10 -	( #	cells	seeded	at	beginning	of	passage	)

Acknowledgements

We thank Professor Pidder Jansen-Dürr (Austrian Academy of Sciences, Innsbruck, Austria) for the HUVEC NS2 cells. We appreciate the technical support of Helle Jacobsen. Financial support was obtained through the FP6 EU project, PROTEOMAGE (LSHM- CT-2005-518230).

A C C E P T E D

M A N U S C R I P T


MD Pascal Lim
email: lim.pascal.hmn@gmail.com

PhD Erwan Donal

PhD Stéphane Lafitte

MD, PhD Genevieve Derumeaux

PhD Gilbert Habib

MD, PhD Patricia Réant

MD, PhD Sophie Thivolet

MD Nicolas Lellouche

MD Richard A Grimm

DO; Pascal Gueret

Multicentre study using strain delay index for predicting response to cardiac resynchronization

Keywords: Cardiac resynchronization therapy, strain, speckle tracking, dyssynchrony

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

INTRODUCTION

Randomized studies have demonstrated that cardiac resynchronization therapy (CRT) has a beneficial impact on reverse remodelling and survival in heart failure (HF) patients with a large QRS interval [START_REF] Abraham | Cardiac resynchronization in chronic heart failure[END_REF][START_REF] Cleland | The effect of cardiac resynchronization on morbidity and mortality in heart failure[END_REF][START_REF] Bristow | Cardiac-resynchronization therapy with or without an implantable defibrillator in advanced chronic heart failure[END_REF][START_REF] Moss | Cardiac-resynchronization therapy for the prevention of heart-failure events[END_REF].

However, up to 40% of patients fail to respond to CRT despite a prolonged QRS duration [START_REF] Cleland | The effect of cardiac resynchronization on morbidity and mortality in heart failure[END_REF][START_REF] Chung | Results of the Predictors of Response to CRT (PROSPECT) trial[END_REF].

Observational studies have consistently shown that the main predictor of responsiveness to CRT is mechanical rather than electrical dyssynchrony. Measurement of regional myocardial electrical-mechanical events using velocity data acquired with tissue Doppler imaging (TDI) has been proposed to enhance identification of mechanical dyssynchrony, and to select patients with wide and narrow QRS duration who may better respond to CRT [START_REF] Bax | Cardiac resynchronization therapy: Part 2--issues during and after device implantation and unresolved questions[END_REF][START_REF] Bax | Cardiac resynchronization therapy: Part 1--issues before device implantation[END_REF]. However, myocardial dyssynchrony markers are usually based on timedelay measurements which are inherently limited, as residual myocardial contraction is not taken into account [START_REF] Lim | Is postsystolic shortening a marker of viability in chronic left ventricular ischemic dysfunction? Comparison with late enhancement contrast magnetic resonance imaging[END_REF][START_REF] Lim | Importance of contractile reserve for CRT[END_REF][START_REF] Foley | Radial dyssynchrony assessed by cardiovascular magnetic resonance in relation to left ventricular function, myocardial scarring and QRS duration in patients with heart failure[END_REF]. Two recent studies demonstrated the limited accuracy of time delay markers to identify responders to CRT in wide (PROSPECT study [START_REF] Chung | Results of the Predictors of Response to CRT (PROSPECT) trial[END_REF]) and narrow (RETHINQ study [START_REF] Beshai | Cardiac-resynchronization therapy in heart failure with narrow QRS complexes[END_REF]) QRS populations.

To overcome these limitations, we recently proposed a new method, the strain delay index [START_REF] Lim | Longitudinal strain delay index by speckle tracking imaging: a new marker of response to cardiac resynchronization therapy[END_REF] (SDI) to predict response to CRT by directly assessing the potential for incremental contractility gain after resynchronization rather than by simply quantifying LV dyssynchrony by regional timing. In the present multicenter study, the SDI was assessed in a large cohort of HF patients referred for CRT (that includes wide and narrow QRS duration) and used to predict future reverse remodelling after CRT.

METHODS

Patients

This study was designed as a multicenter observational clinical trial, and aimed to assess the accuracy of the SDI in identifying response to CRT. Four French university centres equipped with GE ultrasound echocardiography systems participated in the study. To test our hypothesis, we enrolled optimally treated heart failure patients who had undergone implantation of a CRT device. Before CRT, all patients underwent a comprehensive echocardiography examination using a GE ultrasound system. Patients with chronic atrial fibrillation were not included. Before enrolment, all patients gave informed consent and our local ethics committee approved the study. Overall, 235 heart failure patients (120 from Creteil, 65 from Bordeaux-Marseilles, 40 from Rennes and 10 from Lyon) who received CRT were enrolled between November 2008 and November 2009. However, 20% of patients were excluded because a limited acoustic window prevented 2D strain analysis (>2 non analyzable segments). Of the 189 analyzable patients (Table I), 51 (27%) had QRS130ms (narrow) and 138 had a wide QRS complex. In the wide QRS group, all patients had severe LV dysfunction (LVEF ≤35% by 2D echocardiography, mean 25±8%) and had remained symptomatic despite optimal medical treatment (New York Heart Association, NHYA class>II).

Medical treatment included beta-blockers (79%), angiotensin inhibitors (87%) and aldosterone antagonists (40%). In the narrow QRS population, the decision to add biventricular pacing to a standard implantable cardioverter defibrillator was taken by the referring cardiologist, following identification of significant mechanical dyssynchrony. For all patients, the medical decision was independent of the strain delay value.

2D echocardiography:

A comprehensive echocardiographic study (GE, Vingmed System 7, Horten, Norway) was performed before and six months after implantation of the CRT device. During the echocardiographic study, apical views (2-3-, and 4-chamber) with high frame rates (≥50 frames/s) were acquired during breath hold and stored in cine-loop format. All echocardiography data were collected and analyzed in a centralized core lab (Henri Mondor Hospital). The LV volumes and ejection fraction were calculated from the apical 2-and 4chamber views using Simpson's rule. Speckle tracking analysis was performed by an independent operator unaware of the clinical and echocardiographic data.

Strain delay index and LV dyssynchrony:

The method to compute SDI has been described in detail previously [START_REF] Lim | Longitudinal strain delay index by speckle tracking imaging: a new marker of response to cardiac resynchronization therapy[END_REF]. In dyssynchronous ventricles, delayed segments do not contribute fully to end-systolic function. The wasted energy per segment caused by dyssynchrony can be expressed mathematically as the difference between peak (ε peak ) and ES strain (ε ES ). Theoretically, this difference (ε ES -ε peak ) increases with the severity of dyssynchrony (Figure 1). The wasted energy is expected to be greater in the segment with preserved contractility than in the segment with minimal or no residual contractility at a similar degree of delayed contraction. The different steps of SDI processing can be summarized as follows: first, a global strain (ε) curve representing LV function is obtained by averaging 16 regional LV strain curves. Next, the time to the peak of this global-ε curve is used to determine the timing of end-systole (ES) and to compute the strain in ES in the 16 segments. Then, peak and time to peak-ε in the 16 segments is defined as the minimum strain value during the cardiac cycle.

Finally, the difference (ε peak -ε ES ) in each segment (all 16 segments) is summed to generate the SDI (Figure 2). For segments that exhibit positive strain or biphasic strain with a peak positive strain greater than the maximum absolute negative strain, the term (ε peak -ε ES ) is entered as zero for the calculation of SDI. This methodology is based on previous data demonstrating that dyskinetic segments with a predominant stretching motion are unlikely to contribute to CRT response [START_REF] Carasso | Left ventricular strain patterns in dilated cardiomyopathy predict response to cardiac resynchronization therapy: timing is not everything[END_REF]. Time to peak longitudinal strain by speckle tracking is used to calculate the 12-segment (base and mid) SD of time to peak-ε (12SDε). In segments with positive or biphasic strain curves, time to minimum ε was chosen to compute 12SD-ε.

The entire process was computed automatically on an Excel file.

Device implantation:

Implantation of the CRT device was performed in the standard fashion with 3 trans-venous leads inserted.

The right atrial and ventricular (apical site) leads were positioned conventionally. The LV lead was inserted through the coronary sinus and positioned into the lateral or posterolateral cardiac vein. Epicardial implantation (n=3) was required when coronary sinus catheterization failed. Biventricular pacing devices used included those manufactured by Medtronic (Minneapolis, MN, USA), St Jude Medical (Sylmar, CA, USA), and Guidant-Johnson & Johnson (Boston, MA, USA After implantation, the atrioventricular interval was adjusted for optimal diastolic filling (by pulsed mitral Doppler): no VV delay setting was performed. In all patients, a device control was performed at 6 months to confirm the appropriate biventricular pacing (pacing time >80% in all patients included).

Follow-Up:

Clinical characteristics at baseline and 6-months follow-up were obtained from medical reports. The echocardiographic examination was repeated 6 months after CRT. Response to CRT was defined by significant LV reverse remodelling (LVES volume reduction [ESVR] >15%) 6 months after CRT. Changes in SDI and mechanical dyssynchrony after CRT were available in 82 patients that underwent the 6-month follow-up echocardiography with the GE system.

Statistical analysis:

All continuous variables are presented as meanSD and dichotomous data as percentages. To compare numerical data between two groups paired and unpaired Student t-tests were used as appropriate.

Correlations between variables were assessed using a Pearson's linear correlation. Inter and intra-observer reproducibility was expressed as the difference between 2 measurements divided by the mean of these measurements. Sensitivity, specificity and positive and negative predictive values of the SDI were computed using a previously validated cut-off value (SDI >25%). To identify predictors of CRT response, multivariate analysis using the logistic regression model was performed with the inclusion of variables with p value <0.1 by univariate analysis. Statistical significance was considered when the p value was <0.05.

RESULTS

Characteristics of the whole population are summarized in the Table 1. During the 6 month follow up, 65% of patients were clinically improved (change in NYHA score ≥I), while 35% remained stable or worsened their symptoms (3 died, 5 underwent heart transplantation, 5 were hospitalized for recurrent HF and 4 had severe ventricular arrhythmia). The increase in LVEF averaged 11±20% and mean ESVR was 18±25%.

Response to CRT (ESVR>15%) was observed in 114 (60%) patients, and was less in ischaemic (44% vs. 68%, p=0.001, Table 2) and narrow QRS patients (49% vs. 65%, p=0.04, Table 3). Responders to CRT had more preserved global strain (-8±3% vs. -7±3%, p=0.01), more dyssynchrony (12SD-ε=112±37ms vs. 97±30ms, p=0.005) and had greater SDI (39±10% vs. 24±9%, p<0.0001, Table 1). A poor correlation with reverse remodelling after CRT was observed for baseline QRS duration (r=0.17 p=0.02), 12SD-ε (r=0.18, p=0.01) and global strain (r=0.20, p=0.005). In contrast, SDI correlated with LVEF improvement (r=0.45, p<0.0001), ESVR (r=0.61, p<0.0001, Figure 3) and was greater in patients who were clinically improved (36±12% vs. 30±12%, p=0.01). Multivariate analysis adjusted to QRS duration, LV dyssynchrony, global strain, and cardiac aetiology showed that only SDI (OR=1.8, p<0.001) was predictive for response to CRT.

Using the cut-off value previously reported, SDI>25% identified 92% of responders with a positive predictive value of 80% and a negative predictive value of 84% (Figure 4).

Strain delay index in ischaemic patients:

Compared to the non-ischaemic group, ischaemic patients responded less to CRT despite a similar level of dyssynchrony and QRS duration (Table 2). However, compared to the non-ischaemic group, SDI was lower in ischaemic patients (30±11% vs. 35±12%, p=0.01) and differed between responders and nonresponders (Figure 5), while LV dyssynchrony failed to predict response to CRT. Importantly, SDI similarly correlated with ESVR in both the ischaemic (r=0.55, p<0.001) and non-ischaemic group (r=0.61, p<0.0001, Figure 3).

Strain delay index in narrow QRS patients:

Despite a similar level of LVEF impairment, narrow QRS patients in comparison to wide QRS patients had a less dilated ventricle, less severe myocardial dyssynchrony (97±30ms vs. 110±37, p=0.02) and lower SDI value (30±10% vs. 34±13%, p=0.01) (Table 3). Response to CRT was observed in 49% of patients in the narrow QRS group vs. 65% in the wide QRS group (p=0.04). Response to CRT was not associated with LV dyssynchrony (Figure 5) and baseline follow-up QRS duration. In contrast, SDI was greater in responders and correlated with ESVR similarly in the narrow and wide QRS populations (Figure 3).

Changes in dyssynchrony and strain delay index:

In responders, decrease in LV dyssynchrony (-13±48ms vs. +9±49ms, p=0.05) and in SDI (-7±18% vs. 4±11%, p=0.002) was more marked than in non responders. Changes in SDI correlated with ESVR (r=0.35, p=0.001) and LVEF improvement after CRT (r=0.27, p=0.01). In non responders, LV dyssynchrony (+9±49ms) and SDI (4±11%) remained unchanged despite efficient pacing.

Reproducibility:

Intra and inter-observer variability performed in 10 random subjects was 8% and 12%, respectively for automatic 12SD-ε and 5% and 6%, respectively for automatic SDI.

DISCUSSION

Reverse remodelling of the LV after CRT provides an objective measurement of CRT response, that correlates with improved survival and LV function. In this multicenter non randomized study which included a large number of patients, CRT response (ESVR>15%) was observed in 65% of patients with prolonged QRS duration, and in 49% of patients with narrow QRS duration (p<0.05). Left ventricular mechanical dyssynchrony (12SD-ε) derived from time-delay measurement with speckle tracking analysis was lower in non-responder patients but its accuracy is limited to identifying responders to CRT. In contrast, SDI, which quantifies the amount of wasted energy due to LV dyssynchrony, was found to correlate closely with LVEF improvement and reverse remodelling after CRT in both narrow (r=0.59, p<0.001) and wide QRS patients (r=0.66, p<0.0001, Figure 3).

Most CRT studies have demonstrated that up to 40% of patients are non-responders, presumably because mechanical dyssynchrony was absent despite a wide QRS duration [START_REF] Abraham | Cardiac resynchronization in chronic heart failure[END_REF][START_REF] Chung | Results of the Predictors of Response to CRT (PROSPECT) trial[END_REF]. Several methods based on the time-delay measurement of regional wall motion have been proposed to quantify LV dyssynchrony.

Single-centre studies have demonstrated that the 12SD of peak velocity by TDI ( 14) (SD-TDI) and opposing-wall delay ( 15) by radial strain or longitudinal velocity, correlate to the clinical response to CRT. However, in the PROSPECT ( 5) study, neither the SD-TDI (AUC=0.55, sensitivity=0.77, specificity=0.31) nor the opposing-wall delay by TDI (AUC=0.61, sensitivity=53%, specificity=69%) was predictive of response to CRT. Current updated European Society of Cardiology guidelines do not recommend the use of echocardiographic markers of dyssynchrony for selecting patients for CRT [START_REF] Dickstein | focused update of ESC Guidelines on device therapy in heart failure: an update of the 2008 ESC Guidelines for the diagnosis and treatment of acute and chronic heart failure and the 2007 ESC Guidelines for cardiac and resynchronization therapy. Developed with the special contribution of the Heart Failure Association and the European Heart Rhythm Association[END_REF]. The suboptimal accuracy of echocardiographic markers to predict response to CRT has been attributed to the limited accuracy and reproducibility of the Doppler technique in HF patients [START_REF] Mullens | Tissue Doppler imaging in the estimation of intracardiac filling pressure in decompensated patients with advanced systolic heart failure[END_REF][START_REF] Lim | Impact of left ventricular size on tissue Doppler and longitudinal strain by speckle tracking for assessing wall motion and mechanical dyssynchrony in candidates for cardiac resynchronization therapy[END_REF] and/or the use of a single marker approach.

To overcome these limitations, Lafitte et al [START_REF] Lafitte | Validation of an echocardiographic multiparametric strategy to increase responders patients after cardiac resynchronization: a multicentre study[END_REF] proposed a multi-parametric approach using conventional echocardiographic markers to improve identification of responders. In the same way, several authors [START_REF] Delgado | Assessment of left ventricular dyssynchrony by speckle tracking strain imaging comparison between longitudinal, circumferential, and radial strain in cardiac resynchronization therapy[END_REF] have suggested that it may be possible to better predict response to CRT by quantifying myocardial dyssynchrony with the use of speckle tracking [START_REF] Gorcsan | Relationship of echocardiographic dyssynchrony to long-term survival after cardiac resynchronization therapy[END_REF][START_REF] Tanaka | Dyssynchrony by speckle-tracking echocardiography and response to cardiac resynchronization therapy: results of the Speckle Tracking and Resynchronization (STAR) study[END_REF]) that provides a more accurate assessment of myocardial deformation. However, neither the predictive value of mechanical dyssynchrony derived from speckle tracking [START_REF] Miyazaki | Dyssynchrony indices to predict response to cardiac resynchronization therapy: a comprehensive prospective single-center study[END_REF] for responsiveness to CRT nor the multi-parametric [START_REF] Bordachar | Low value of simple echocardiographic indices of ventricular dyssynchrony in predicting the response to cardiac resynchronization therapy[END_REF] approach have been confirmed. The controversial accuracy of time delay measurement using either TDI or speckle tracking may be explained by the fact that delayed contraction is a nonspecific marker of myocardial dysfunction that can be seen in scar [START_REF] Lim | Is postsystolic shortening a marker of viability in chronic left ventricular ischemic dysfunction? Comparison with late enhancement contrast magnetic resonance imaging[END_REF][START_REF] Foley | Radial dyssynchrony assessed by cardiovascular magnetic resonance in relation to left ventricular function, myocardial scarring and QRS duration in patients with heart failure[END_REF], fibrosis [START_REF] Kawara | Activation delay after premature stimulation in chronically diseased human myocardium relates to the architecture of interstitial fibrosis[END_REF][START_REF] De Bakker | Fractionated electrograms in dilated cardiomyopathy: origin and relation to abnormal conduction[END_REF], and viable ( 27) myocardium and may be further impacted by loading conditions [START_REF] Skulstad | Postsystolic shortening in ischemic myocardium: active contraction or passive recoil[END_REF][START_REF] Park | Impact of loading condition on the 2D speckle tracking-derived left ventricular dyssynchrony index in nonischemic dilated cardiomyopathy[END_REF]. This limitation particularly affects ischaemic patients but may be overcome by investigating myocardial viability and contractile reserve [START_REF] Lim | Importance of contractile reserve for CRT[END_REF][START_REF] Costa | Prospective validation of stress echocardiography as an identifier of cardiac resynchronization therapy responders[END_REF][START_REF] Bleeker | Effect of posterolateral scar tissue on clinical and echocardiographic improvement after cardiac resynchronization therapy[END_REF][START_REF] Hummel | Extent of myocardial viability predicts response to biventricular pacing in ischemic cardiomyopathy[END_REF] as a complement to LV dyssynchrony to better identify responders.

In a recent study [START_REF] Lim | Longitudinal strain delay index by speckle tracking imaging: a new marker of response to cardiac resynchronization therapy[END_REF], we proposed a new method, the SDI that combines within a single marker mechanical dyssynchrony measurement and myocardial contractility, to predict response to CRT. Strain delay index is not a simple measurement of contractility or time delay but a combination (and relative weighting) of both of these parameters (Figure 1). The importance of considering contractility and dyssynchrony to predict responsiveness to CRT is particularly crucial in ischaemic patients [START_REF] Carasso | Left ventricular strain patterns in dilated cardiomyopathy predict response to cardiac resynchronization therapy: timing is not everything[END_REF][START_REF] Reant | Identification and characterization of superresponders after cardiac resynchronization therapy[END_REF]. Indeed, ischaemic patients respond less to CRT despite a similar level of dyssynchrony and QRS duration than non-ischaemic patients. This may be explained by the presence of extensive scar in delayed segments that prevents efficient biventricular pacing [START_REF] Lim | Importance of contractile reserve for CRT[END_REF][START_REF] Costa | Prospective validation of stress echocardiography as an identifier of cardiac resynchronization therapy responders[END_REF][START_REF] Bleeker | Effect of posterolateral scar tissue on clinical and echocardiographic improvement after cardiac resynchronization therapy[END_REF][START_REF] Hummel | Extent of myocardial viability predicts response to biventricular pacing in ischemic cardiomyopathy[END_REF]. In ischaemic and non ischaemic patients, leads positioned in scar [START_REF] Bleeker | Effect of posterolateral scar tissue on clinical and echocardiographic improvement after cardiac resynchronization therapy[END_REF] or non desynchronized segments [START_REF] Ypenburg | Optimal left ventricular lead position predicts reverse remodeling and survival after cardiac resynchronization therapy[END_REF] and severe myocardial dysfunction related to myocardial fibrosis and remodelling may also limit the beneficial effects of CRT. Reant et al [START_REF] Reant | Identification and characterization of superresponders after cardiac resynchronization therapy[END_REF] recently demonstrated that responders and super responders to CRT had greater longitudinal global strain than non responders. Scar and fibrosis in delayed segments may be investigated by stress echocardiography or magnetic resonance and nuclear imaging. However, the physiological approach of SDI may provide a simple alternative to these complex and costly investigations. The SDI incorporates the degree of impaired contractility related to myocardial dyssynchrony rather than dyssynchrony alone. This explains the close correlation between strain delay index and reverse remodelling after CRT in ischaemic and non-ischaemic patients. In the present study, we demonstrated that SDI and decrease in SDI after CRT was greater in responders. In addition, decrease in SDI after CRT appears to be correlated with reverse remodelling and LVEF improvement. These results suggest that the SDI may be inferred as the gain of contractility expected after optimal resynchronization therapy.

Narrow QRS patients may respond to CRT when mechanical dyssynchrony exists [START_REF] Van Bommel | Association of intraventricular mechanical dyssynchrony with response to cardiac resynchronization therapy in heart failure patients with a narrow QRS complex[END_REF]. In the present study, we found that narrow QRS patients who responded to CRT had a similar SDI to wide QRS patients despite lower mechanical dyssynchrony. This can be explained by the fact that the narrow QRS patients had less advanced heart disease and more preserved contractility in delayed segments. Indeed, the wasted energy in segments with limited dyssynchrony (5-10% of RR interval) and preserved contractility can reach a similar level to severely desynchronized segments (>10% delay from ES, Figure 6). Finally, these results highlight the concept of mechanical dyssynchrony by demonstrating that limitation of time delay markers is not simply related to the accuracy of the method used to quantify timing but rather to the need to consider contractility. The SDI may overcome this limitation by quantifying the wasted energy related to dyssynchrony that can be recruited by biventricular pacing. The wasted energy is expected to be higher in less advanced heart failure patients (NYHA I-II or/and LVEF>35%) with wide QRS complex. In this population, recent studies have demonstrated a substantial benefit of CRT [START_REF] Lubitz | Effectiveness of cardiac resynchronization therapy in mild congestive heart failure: systematic review and meta-analysis of randomized trials[END_REF][START_REF] Chung | Cardiac resynchronization therapy may benefit patients with left ventricular ejection fraction >35%: a PROSPECT trial substudy[END_REF]. The use of SDI may be particularly useful for selecting patients who will most benefit from CRT, especially when the QRS complex is moderately enlarged (<150ms).

Study Limitations: Despite the large number of patients included in the study, only a limited number of patients with a narrow QRS duration underwent CRT. The size of the study is a limitation to firmly concluding that the SDI in narrow QRS patients is an accurate method for identifying responders to CRT.

A randomized study should be performed in this setting to confirm our results. In addition, the follow up period should be extended to fully understand reverse remodelling related to CRT and its association with SDI. Further studies should be performed to assess the impact of myocardial scar extent and lead position relative to scar, on the accuracy of SDI in predicting response to CRT. Radial and circumferential strains were not used to quantify SDI because their values cannot be computed from the whole myocardium and reproducibility remains limited. Finally, the study excluded 20% of patients with limited acoustic windows to provide the real accuracy of SDI.

Conclusions:

The present multicenter prospective study, that included non-selected patients referred for CRT suggests that SDI, which considers LV dyssynchrony and residual contractility, may identify responders to CRT.

Figure 1 :

 1 Figure 1: Changes in wasted energy (ε peak -ε ES ) according to myocardial dyssynchrony and contractility.

Figure 2 :

 2 Figure 2: Strain delay index is defined as the sum of the wasted energy ( ESpeak ) caused by LV

Figure 3 :

 3 Figure 3: Correlation between strain delay index and ESVR according to cardiac aetiology (A) and to QRS

Figure 4 :

 4 Figure 4: ROC curves showing the accuracy of strain delay index in predicting response to CRT. SDI =

Figure 5 :

 5 Figure 5: Left ventricular dyssynchrony was greater in responders with a wide QRS duration and non

Figure 6 :

 6 Figure 6: Data from wide (A) and narrow (B) QRS patients demonstrating that the wasted energy per

Sources of Funding:

This work was funded by the French Federation of Cardiology and the French Society of Cardiology.

Conflict of Interest Disclosures: none


Konrad Hinsen
email: konrad.hinsen@cnrs-orleans.fr

A data and code model for reproducible research and executable papers

published or not. The documents may come

Introduction

One of the hallmarks of science is reproducibility: a scientific study must be documented to the point that another researcher can follow the same steps and obtain the same results. Computational science currently falls short of this goal because the programs and the input data sets for a computation are rarely 1 published. Recent efforts to improve this situation by creating suitable tools and awareness of the problem in the scientific community have been conducted under the somewhat provocative label "reproducible research" [START_REF] Claerbout | Guest editors' introduction: Reproducible research[END_REF]. Executable papers carry this idea one step further: a single published digital object should contain an article describing a scientific study and its results as well as all the data and program code necessary for repeating the calculations and visualizing the results. Furthermore, this digital object must be suitable for long-time archival and for inspection on a large variety of computing platforms, present and future.

The work described here concentrates on the infrastructure required for making computational research reproducible and publishable. It addresses the question of how data, code, and human-readable text can be combined and stored efficiently in a single file that can be used on current and future computer architectures. The proposal uses already available technology as much as possible, pointing out where shortcomings of these technologies are likely to require improvements in order to make the proposed executable paper format suitable for production use in a wide range of computational science domains. Key features of the proposed format are

• platform independence through virtual machine technology

• automatic verification of all computational results for which code is provided by the authors, and authentication by electronic signatures of the results for which no code is provided

• protection against viruses and other malicious code

• provenance tracking for computational and experimental data with optional electronic signatures

• references to items in other executable papers allow re-use of previously published data and code An aspect not addressed by this proposal is user interfaces, because it is important to clearly separate data representation and user interfaces. User interface choices depend on both personal preferences (e.g. some people prefer graphical user interfaces, others command line interfaces) and the hardware being used (a desktop computer and a smartphone have very different requirements). Data representation, on the contrary, must be as uniform and precisely defined as possible to guarantee portability and longevity of electronic documents.

A key design principle behind this proposal is the notion that "code is data". There is in fact no fundamental difference between executable programs and other data items: they are stored in memory or in files and can be transformed or analyzed by other programs. For both code and other data, clarity, portability and long-term usability require well-defined and carefully designed data models. The core of this proposal is exactly that: a data model for reproducible research and for executable papers. The next section discusses the most important choices to be made: how to represent code and other data, and how to package data and metadata to make an executable paper.

2 Representation of data, code, and text

File storage

For reasons of integrity and for simplifying automated treatment, an executable paper must consist of a single file. The file format chosen must thus be able to accommodate raw data (including potentially very large data sets), metadata (for provenance tracking and for use in data models), executable code, text including markup for equations etc., and machine-treatable references to other executable papers (for corrections, updates, comments, etc.).

The Hierarchical Data Format 5 (HDF5) [START_REF] The | Hierarchical data format version 5[END_REF] fulfills all these requirements. It is a format specifically designed for scientific computing and addresses issues specific to this domain, e.g. the prevalence of large homogeneous data sets such as arrays. HDF5 is already widely used for implementing data models for experimental data (e.g. NeXus [START_REF] Klosowski | NeXus: A common format for the exchange of neutron and synchrotron data[END_REF], used in neutron and X-ray scattering) or computational results (e.g. CGNS [START_REF] Poirier | The CGNS system[END_REF][START_REF]CFD general notation system (CGNS[END_REF], used in Computational Fluid Dynamics).

An HDF5 file has an internal tree structure that resembles a file system as used by conventional operating systems, with group nodes playing the role of directories and datasets (leaf nodes) storing the actual data. Data sets can be simple values but more typically are arrays or tables of values. Data sets can be of fixed or variable size and use optional features such as compression. Individual values are defined by their data type, which can be a built-in standard type such as "16-bit unsigned integer" or "variable-length character string" or a userdefined (compound) type. There is also an "opaque" type for storing arbitrary byte sequences which is suitable for storing binary code representations. The wide range of data types and dataset types provides the flexibility required for storing any kind of data and choose a suitable representation that is the right compromise between file size, access time, and convenience.

Each group and dataset in an HDF5 file can have metadata attached in the form of named attributes. Attribute values can be of any HDF5 data type. A data type of specific interest for attribute values is the reference type, which stores a pointer to a dataset. Reference-type attributes can be used to store dependency information between datasets.

The HDF5 format is implemented as a highly portable library written in the C language. It is developed and maintained by the HDF Group, located at the University of Illinois. Interface layers for C++, Fortran, and Java (using the Java Native Interface (JNI)) are provided by the HDF group as well. Interface layers for a large number of other programming languages are available from other sources. A large number of generic tools for transforming and visualizing data stored in HDF5 files is available from the HDF Group and others. Two important more recent developments are worth mentioning here because they ensure the utility of HDF5-based formats on future computing platforms. A parallel-I/O version of HDF5, available from the HDF Group, allows its efficient use in supercomputing environments. The OPeNDAP software [START_REF]OPeNDAP: Open-source Project for a Network Data Access Protocol[END_REF] provides remote access to datasets stored in HDF5 (and other) files. Remote access is important for handling very large datasets that cannot be transferred completely to every user's computer, but also for browsing tools on portable devices such as smartphones and tablet computers.

Domain-specific data models

Many kinds of scientific data require non-trivial and non-obvious representations. For example, storing the configuration of a protein implies storing the positions of all atoms as well as sufficient information about the chemical structure of the protein to be able to identify the structural role of each atom as well as its most important chemical properties. The correct interpretation of such data relies on a well-defined data model.

Data models are necessarily domain-specific, and therefore cannot be part of a universal format for executable papers, with the exception of simple and frequently used data models e.g. for time series or mathematical equations. Each dataset or group of datasets has a metadata tag (an HDF5 attribute whose name is fixed by convention) indicating the data model to be applied.

It is expected that each electronic journal will select a set of data models relevant for its community and apply verification tools for those data models upon submission of an electronic paper. It is also expected that domain-specific support libraries will emerge that facilitate the use of HDF5-based data models and spare scientists and scientific software developers the effort of using the rather complex HDF5 programming API directly.

Dummy data sets

If a dataset can be reconstructed from other information in an executable paper, why store it at all? The main reasons are

• Time: recalculation of the data may be feasible, but too slow for interactive browsing.

• Space: the data set may be too big to be stored in memory, requiring on-disk storage prior to visualization.

• Simplicity of access: analysis or visualization software may not have the functionality to recompute the data.

In these situations, the explicit presence of a dataset is a form of caching. For long-term storage, where data size is the most prominent criterion, it may be preferable to delete such data sets and replace them by dummy data sets that retain only the metadata, permitting recalculation at a later time.

Executable code

The representation of code inside an executable paper is the most difficult aspect because of conflicting requirements. A suitable code representation must be platform-independent and stable over very long time spans in order to permit its use on future computer systems whose characteristics are still unknown. It should rely on a run-time system that can block erroneous or malicious behavior. It should have an unambiguous specification of the outcome of arbitrary computations. It should admit a wide range of programming languages, current and future. It should not stand in the way of efficient program execution on a wide range of computer architectures. Finally, it should permit the re-use of existing scientific software as much as possible. In the following I will examine existing code representations according to these criteria.

Source code

The availability of source code for computational procedures in an executable paper is highly desirable because it is the most precise documentation of the algorithms, and it permits readers to explore modifications. However, source code is not a good choice as the primary code representation in an executable paper. First of all, there is a very large number of programming languages in use, of which many have machine-or vendor-specific dialects. An executable paper format based on source code would have to allow only a well-defined (and thus fixed) set of programming languages, which would be an obstacle to future progress in software development. Moreover, many programming languages are intentionally ambiguous, leaving the precise meaning of some constructs undefined in order to allow compilers to apply machine-specific optimizations.

Floating-point arithmetic is particularly concerned by this problem. As a consequence, the same source code compiled on different machines can yield different results, which violates the reproducibility requirement. Finally, the use of source code would require a runtime system for working with executable papers to include compilers or interpreters for all allowed languages.

Machine code

Machine code has the advantage of being unambiguous and able to exploit all features of a computer at maximum performance.. However, it is also inherently unportable, being specific to a processor or a processor family and in practice also to an operating system or even a precise version of an operating system. Machine code also poses security problems as it is very difficult to prevent malicious behavior. Some of these problems can be circumvented by making a complete software installation, including the operating system, available through virtualization technology and run it in a secure environment. Moreover, emulators can be used to execute machine code for one machine on a different machine. However, both approaches lead to performance loss (which is severe in the case of emulators) and are not suited for long-term archival because computer architectures will inevitably change in the future.

Virtual machine bytecode

Virtual machines executing bytecode (also called p-code) provide an intermediate representation between source code and machine code that is both portable and unambiguous. They also make it easy to restrict access to resources such as local files for security reasons. The main disadvantage of using a virtual machine is the loss of performance compared to optimized machine code. However, Just-In-Time (JIT) compilers have reduced the performance gap between virtual machines and native code, and this trend is likely to continue.

There are at the moment two widely used virtual machines that have good enough performance and sufficient programming language support to be candidates for use in an executable paper format: the Java Virtual Machine (JVM) [START_REF] Lindholm | The Java Virtual Machine Specification[END_REF] and the Common Language Infrastructure (CLI) [START_REF]ECMA Standard 335: Common Language Infrastructure CLI[END_REF]. In the following I will concentrate on the JVM because I am more familiar with it. However, the CLI might well be an equally good or even better choice.

The JVM is not an optimal choice for scientific computing. It wasn't designed for this field and therefore has a couple of shortcomings. It lacks value types (types whose data fields are inlined into a parent object or an array and don't require a separate object allocation with the associated garbage collection overhead), which would make many scientific data items (complex numbers, geometric data such as points or triangles, ...) more efficient in terms of memory use and CPU time. The JVM's handling of floating-point arithmetic has been criticized frequently because of its lack of some important IEEE 754 features. In modern JVMs, floating-point arithmetic is not platform-independent unless "strict math" is explicitly requested by the program; this feature was introduced mainly for improving performance on the popular Intel x86 architecture.

Perhaps the most significant obstacle to rapid adoption of any JVM-based infrastructure for executable papers is the lack of production-quality compilers for the most popular scientific programming languages (Fortran, C, C++) that produce JVM bytecode. However, proof-of-concept compilers for C [START_REF] Given | Clue: an ANSI C compiler targeting high level languages[END_REF] and Fortran [START_REF] Fox | A prototype of Fortran-to-Java converter[END_REF][START_REF]f2j, a Fortran to Java compiler[END_REF] already exist. The use of an x86 emulator such as JPC [START_REF] Preston | Pure Java x86 emulator[END_REF] may be of interest as a temporary solution for running legacy software, in spite of the significant performance cost.

Considering the numerous advantages of virtual machines (portable binary representation, portable mixed-language programming, automatic memory management, code optimization based on run-time information) and the fact that modern JIT compilers have almost eliminated their performance overhead, I expect the scientific computing community to adopt virtual machines in the long run for most computational tasks. Executable papers may well be one important motivation to move in this direction. New compilation tools such as LLVM [START_REF] Lattner | LLVM: A compilation framework for lifelong program analysis & transformation[END_REF] and VMKit [START_REF] Geoffray | VMKit: a substrate for managed runtime environments[END_REF] are an important technological innovation for such a move. They blur the frontier between native machine code and virtual machine bytecode, allow a better integration between the two worlds, and facilitate the design and implementation of virtual machines. In the not-too-distant future, scientists will generate portable bytecode for publication and highly optimized native code for their supercomputer from the same source code.

Source code compilation

Even with virtual machine bytecode as the primary code representation in an executable paper, it is highly desirable for its users to have access to the source code as well, because it is the only complete description of the algorithms being used. This raises the question of ensuring that source code and bytecode are equivalent such that readers can trust that the source code they read corresponds to the computations they run.

The problem is exactly the same one as for dependencies between data sets, and it can be solved using the same approach: bytecode datasets are tagged with a reference to the source code dataset from which they can be recovered by compilation. This requires only that compilers be integrated into the executable paper ecosystem in the same way as computational software is, i.e. compilers (and associated tools) must be available in the form of virtual machine bytecode. A few JVM-based compilers already fulfill this requirement (see e.g. [START_REF] Unkrig | Janino, a super-small, super-fast Java compiler[END_REF], [START_REF]Groovy, an agile dynamic language for the Java platform[END_REF], [START_REF] Hickey | Clojure, a dynamic programming language that targets the Java Virtual Machine[END_REF]), and I expect their number to grow. There is also a unified compilation API for the JVM [START_REF]Apache commons JCI, a Java compiler interface[END_REF]. Once again, the adoption of a virtual-machine-based executable paper framework is likely to create the necessary motivation for writing compatible compilers.

Scripts

The above discussion has concentrated on the traditional way of writing scientific programs: source code is translated to bytecode or machine code, which is then run. Another mode of operation has been gaining popularity recently: the use of scripting languages at the highest levels of problem specification, making use of libraries written in lower-level languages for efficiency. Examples of popular scripting languages in scientific computing are Matlab, Python, and R. These languages are typically interpreted, with the interpreter proposing an interface to compiled low-level code.

In the proposed framework for executable papers, a scripting language can be made portable and unambiguous by including its interpreter (usually in the form of a reference, see section 2.7) in the paper itself as JVM bytecode. The use of scripting languages should even be particularly encouraged, because it makes the higher-level algorithmic logic of the computation clearer and easier to modify by the reader. Browsers for executable papers should include a script editor that permits readers to explore the influence of parameters and of algorithmic choices.

Text, tables, and figures

A traditional print-based scientific publication has as its prime ingredients text (including mathematical formulas), tables, and figures. Tables and figures are replaced by dynamically generated table views (i.e. code) and visualization scripts in executable papers. As for text, several useful representations are already in use, and tools for working with them are widely available. The transition to executable papers requires only one new feature: the possibility to include a reference to datasets and executable programs. This feature is readily provided by established hyperlink mechanisms.

The choice of a suitable text representation depends mainly on issues outside of the scope of this article: the needs of authoring tools and browsers. An obvious criterion is the suitability of the representation for algorithmic treatment, which excludes purely visual formats such as PDF. Candidates are light-weight content-oriented markup languages such as reStructuredText [19] or Markdown [START_REF]Markdown[END_REF], complemented by an equation language, XML-based formats such as Doc-Book [START_REF]The DocBook schema version 5.0[END_REF], or traditional typesetting languages such as T E X [START_REF] Ervin | The TeXbook[END_REF].

An aspect that deserves a special discussion is the representation of mathematical equations, which are widely used in many branches of science. Often they are used exclusively for the benefit of the human reader, for defining quantities and for documenting computational methods. They then have the same role as text and should be included in the latter using suitable markup. Equations that correspond directly to a program should be accompanied by a reference to the latter. It is, however, unrealistic to try to establish a direct correspondence between a program and an equation that documents it. In most practical situations, a program is described by more than one equation plus a context (definition of the quantities) provided in text form. Moreover, the algorithm implemented in the program usually deviates from the written equations for reasons of efficiency or data representation.

Mathematical equations take a more active role in publications in which their manipulation (in the form of symbolic computation) is an essential part of the paper's content. They should then be considered data and represented by a domain-specific data model just like other data. Their graphical display is a special case of visualization, to be handled by visualization scripts. It should be noted that the World Wide Web Consortium's MathML markup language [START_REF]The World Wide Web Consortium[END_REF] recognizes the two distinct roles of mathematical equations by defining separate markup languages for the display and the semantics of mathematical equations.

Visualization scripts

The most general specification for a visualization (defined here as any dynamically generated data display, i.e. including tables) is a program, meaning that the executable paper specification does not need to provide any special visualization features. However, it would be inconvenient to leave the choice of visualization routines fully to the authors of an executable paper. This would represent a burden for the authors and lead to a lack of coherent user inter-face for readers. Visualization should be handled by the run-time system for working with executable papers (see section 4) and by its domain-specific extensions. Ideally, the run-time system provides a special scripting language for visualization which is then used in papers to specify interactive data displays.

In traditional print publishing, papers cite other papers using human-readable text, usually following some formatting convention. Such citations are notoriously difficult to treat electronically because of format variations. With the advent of electronic publishing, Digital Object Identifiers (DOIs) have replaced text-based references for use in databases and other electronic resources [24]. A DOI is a character string that uniquely identifies an electronic object. An organization issuing a DOI is responsible for maintaining access to the object and its associated metadata. DOIs are not limited to readable documents; for example, the DataCite consortium [START_REF]The DataCite consortium[END_REF] promotes the use of DOIs for access to scientific data. DOIs are clearly a suitable way to refer to executable papers as well.

However, executable papers should permit more fine-grained references than DOIs in order to target specific datasets inside another executable paper. Such a mechanism permits the re-use (with proper attribution) of data and code published earlier. It can also be of interest to publish an executable paper explicitly in several pieces, i.e. to facilitate sharing of very large data sets among several publications. Since inside an HDF5 file, every dataset is uniquely identified by its path (the sequence of groups to be traversed in order to reach the dataset), which can be represented by a character string, a reference into another executable paper consists of two character strings: the DOI of the paper, and the path of the dataset inside the paper.

Provenance tracking

Data sets that are computed via code provided in an executable paper itself need no provenance information, as they can be recomputed at any time for verification. However, every paper must include information that cannot be regenerated automatically, as otherwise the paper would be of no interest. At the very least, a paper contains program source code and presentation text that are original work of its authors. Experimental data and input parameters for computations also fall into this category. Another category of data that cannot be reproduced is data taken over from legacy sources (print publications, databases) that cannot be referenced directly. Finally, there are data items that could in principle be reproduced but with means not easily available to every reader. This category covers computations that require a supercomputer, but also results obtained from proprietary or legacy software that cannot be included in an executable paper.

All the data items cited above should be accompanied by provenance information, and provenance information should be verifiable. A suitable verification scheme for all the cases cited is an electronic signature, i.e. a digital fingerprint of the data encrypted using a public-key infrastructure. In the case of original work, the electronic signature identifies the authors, whereas non-reproducible computational results could be signed automatically by the batch processing system of a supercomputing center, using a signature identifying the computer system being used.

The digital fingerprints from the signatures can also be used to set up a database of published datasets, indexed by fingerprints, and consulted to detect basic forms of plagiarism. However, more sophisticated detection systems for plagiarism will be required, as it is rather straightforward to copy someone else's dataset and modify it slightly such as to change the fingerprint. In the preceding sections, I have explained the important design choices of my proposal for an executable paper format, but less relevant details were intentionally left out. In this section, I briefly present a concrete data layout to convey a clearer impression of what an executable paper looks like. Figure 1 illustrates the dependencies between data items in an executable paper.

Executable papers

An executable paper consists of a single HDF5 file. Published executable papers are identified by a Digital Object Identifier (DOI) and contain their own DOI in their metadata, permitting automated authenticity checks.

The file's root group contains three subgroups "data", "code", and "text". The "data" group contains groups and datasets that each have a metadata tag (HDF5 attribute) specifying the domain-specific data model they respect. Other attributes contain provenance information (for non-recomputable data), or references to dependencies plus a reference to the program (in the "code" subgroup) that re-calculates the dataset from its dependencies. There can also be a "legend" attribute pointing to a plain-text explanation for human readers.

The "code" subgroup contains any number of HDF5 opaque datasets storing a jar file per dataset. Attributes of these datasets store relevant metadata: provenance, reference to the source code, reference to the compilation script run to rebuild the jar file. The "code" subgroup also contains any number of scripts. A script is a character-string dataset whose attributes store references to the jar files making up the interpreter, and the JVM command for running the interpreter. Finally, the "code" subgroup can contain any number of source code files and groups containing source code files. A source code file is stored in a character-string dataset with provenance information in the metadata.

The "text" subgroup contains any number of character-string datasets representing text using a markup language. Each dataset has provenance information in its metadata. One dataset must be called "main" and contains the main presentation text that is shown to the reader upon opening the executable paper. References (hyperlinks) in each text can refer to other text items or to scripts/programs that produce visualizations. References can also point to datasets in the "data" section, allowing the user to open a data browser window for a specific dataset.

It is important to recall that references can point to items in previously published executable papers, using a combination of the DOI and the data item path inside the file. In fact, executable papers cover a much wider range of publishable items than traditional printed papers. Some important special cases of executable papers are:

• Raw data sets, usually from experiment or observation. Such papers contain data and text (for documentation), but no code.

• Program libraries, containing code and text (for documentation) but no data.

• Comments on published papers, containing only text and optionally visualization scripts.

Such items, which currently do not fit into the scientific publication system, become normal publications that have authors, a publication date, and which can be cited. They are automatically integrated into the archival systems operated by scientific publishers, ensuring their long-term conservation in unmodified form, and thus the reproducibility of other papers that refer to them. Publishing such items separately also provides a means of asserting authorship on important contributions even if they are not in themselves scientific results.

Software infrastructure

Working with executable papers requires a variety of software tools. For successful adoption of an executable paper format, it is essential that it can be integrated into as many existing tools as possible, including

• interactive scientific computation environments

• workflow managers

• authoring environments for composing electronic papers

• software development environments

• visualization software

In addition, new functionalities are required that can be either integrated into existing tools or provided by new tools:

• interactive browsers for desktop computers that let readers explore the text and data as well as re-run the computations and modify parameters

• server-side tools with the same functionality, accessed via a standard browser interface by the user working on a desktop computer or a less powerful device such as a smartphone

• automatic verification tools that re-run all computations in the paper and compare the results to the ones supplied by the authors

• authentication tools that verify the electronic signatures

The current proposal is based on two existing and well-established technologies: the HDF5 format and library for data storage, and the Java Virtual Machine for code execution. The HDF5 format is already widely supported by scientific computing tools. The JVM has found its niches in scientific computing, but cannot be considered popular yet. This is partially due to its real or perceived defaults (see section 2.4.3), but mostly due to a lack of scientific libraries and compilers for the most popular scientific programming languages.

The most important support software that needs to be written for using executable papers is a run-time system for executing the code stored in such a paper. This run-time system must obviously include a JVM and the HDF5 library. It must be able to run programs from a paper in a secure environment in which access to local resources on the user's computer is limited to reading data from executable papers and writing data to a single HDF5 file specified by the user. In addition, a good run-time system should provide support libraries for data management, visualization, etc. The generic run-time system would in practice be complemented by domain-specific extensions that handle domain-specific data models and access to domain-specific resources such as public databases.

Related work

The literate programming and reproducible research features [START_REF] Schulte | Active documents with org-mode[END_REF] in Emacs orgmode [START_REF] Dominik | Org-mode for emacs[END_REF] were a major inspiration for this work. They permit the storage of text, code, and data in a single plain-text file and provide an interface between code sections written in different languages. Org-mode is limited, however, by its plain-text representation; data sets cannot become very large and the code contained in such a file necessarily depends on external software (compilers and interpreters, usually also libraries) that cannot be included in the same package.

The SHARE system [START_REF] Van | Supporting the internet-based evaluation of research software with cloud infrastructure[END_REF] also uses virtual machine technology (coupled with network access through a browser) to provide platform-independent and secure access to research code and data. It differs in using system virtualization to store the researcher's entire working environment, including the operating system, the whole file structure, and programs in the form of machine code. The obvious advantage of this approach is the possibility to use existing software and data formats. On the other hand, system-level virtual machines are impractical for long-term storage because of their size and their narrowly defined hardware requirements. Moreover, data and code stored in this way are not reusable in other work.

Various existing scientific computing tools, including workflow managers such as VisTrails [START_REF]VisTrails[END_REF] or Kepler [START_REF]The Kepler Project[END_REF], interactive scientific computation environments such as Matlab [31] or Spyder [START_REF] Raybaut | Scientific PYthon Development EnviRonment (Spyder)[END_REF], or automated logging systems such as Sumatra [START_REF] Davison | Sumatra, a tool for managing and tracking projects based on numerical simulation or analysis[END_REF], handle much of the dependency information that makes the difference between a collection of code and data and an executable paper. Another interesting approach to dependency handling that also permits the unification of code and data is the use of software build tools for computation. An example is the use of SCons [START_REF] Knight | Scons software construction tool[END_REF] in the Madagascar tool suite [START_REF]Madagascar, an open-source software package for multidimensional data analysis and reproducible computational experiments[END_REF]. All of these tools could probably be modified easily to work with the data model described in this work.

Conclusion

The preceding sections propose a data and code representation model that support reproducibility in computational science and the publication and archival of reproducible computational results. It addresses many of the open issues in the design of executable papers: Executability The proposed system permits authors to provide the full code implementing their computational procedures with their publication. It also permits the publication of code libraries that other papers can build on.

Short and long-term compatibility The use of a virtual machine (JVM) ensures code portability on current and future platforms. The HDF5 library ensures the portability of data storage. The long-term compatibility of the format relies on the long-term maintenance of the JVM and the HDF5 library, or on translation mechanisms for converting first-generation executable papers to an eventual future format. Long-term archival by scientific publishers of all data and code guarantees long-term reproducibility of results.

Validation The proposed format permits an automatic validation of all results for which code is provided, and the automatic verification of electronic signatures for all other data.

Universality Domain-specific data models and domain-specific extensions to the run-time system ensure the applicability to all domains of computational science. The use of virtual machine bytecode as the principal code representation admits a wide range of programming languages.

Use of supercomputers Results that cannot be easily reproduced due to exceptional system requirements (i.e. the use of supercomputers) are treated like experimental data: they are supplied with provenance information and an electronic signature certifying their origin. Program code for reproducing them can still be supplied and used in the future if the required computing power becomes more easily available.

Data size The HDF5 library can handle large datasets efficiently and provides support for parallel I/O. Large datasets can be published separately and accessed from analysis papers through references. Remote access to subsets of a large dataset stored on a server eliminates the need to transfer large datasets completely to every user's computer.

Provenance The origin of each data item can be documented by provenance metadata and certified by an electronic signature. Computational steps are documented by their program code and can be reproduced at any time.

Security The use of a virtual machine (JVM) permits the execution of code contained in a paper in a secure run-time environment that protects users against viruses, trojans, and spyware.

The biggest obstacle to the rapid adoption of this proposal is the requirement for all computational code to be based on a virtual machine, which excludes many currently popular programming languages and tools. However, there is no fundamental difficulty with writing the required tools, and the adoption of executable papers by major scientific publishers is likely to encourage such developments.

The implementation of the basic run-time system requires only a modest development effort because most of the technology already exists. Implementing domain-specific extensions represents a much larger effort, because of the large number of domains with specific requirements, and also because many domains of computational science did not yet develop well-defined data models.

Figure 1 :

 1 Figure 1: Data dependencies in an executable paper. The items on the left of the dashed line have no incoming arrows and cannot be recomputed. The script uses the bytecode libraries, but does not depend on them. However, data items that depend on the script indirectly also depend on the bytecode libraries.

Acknowledgements

This work was supported by the Agence Nationale de la Recherche (Contract No. ANR-2010-COSI-001-01).