Emmanuelle Roque D'orbcastel
email: emmanuelle.roque@ifremer.fr

Yvan Bettarel

Marion Dellinger

Bastien Sadoul

Thierry Bouvier

Justin Monin Amandé

Laurent Dagorn

Benjamin Geffroy

Measuring cortisol in fish scales to study stress in wild tropical tuna

Keywords: scales, cortisol, stress, welfare, skipjack tuna, yellowfin tuna

Cortisol is recognized as a physiological indicator of stress in fish. However, this hormone is typically measured in plasma samples. In this study, cortisol content was evaluated for the first time in the scales of tropical tuna. The sample included 20 skipjack tuna (Katsuwonus pelamis) and 25 yellowfin tuna (Thunnus albacares) collected in the Atlantic (the Gulf of Guinea off the Ivory Coast) and the Indian Ocean (off Reunion Island), and their scales were analyzed using an ELISA test. The results indicated that on average, cortisol concentration (1) did not show any geographical pattern, (2) was independent of fish size, and (3) was significantly higher in skipjack (mean: 4.75 ± 6.56 ng g -1 ) than in yellowfin (mean: 1.65 ± 1.85 ng g -1 ), although the difference was mainly due to four skipjack individuals. Larger datasets would be needed to confirm any species difference in cortisol concentration. The particularly high cortisol concentration observed in four individuals may be due to fish from different schools, suggesting the need for future sampling from free-swimming schools to investigate intra-and interschool variability in cortisol concentration. This study opens the door for future research, including the collection of scales on tagged individuals, to investigate the links between chronic stress and behavior in these species.

Introduction

The capacity of fish to cope with environmental challenges is of primary importance in a rapidly changing world. The stress axis plays a central role in the coping capacity of fish, as it diverts energy from growth and reproduction to cope with acute environmental changes [START_REF] Sadoul | 5-Stress and Growth[END_REF].

The two stress-related response mechanisms in fish are brain-sympathetic chromaffin cells (BSC) and the hypothalamus-pituitary-interrenal (HPI) axis. The underlying mechanisms of the stress axis in fish have been well described [START_REF] Barton | Stress in fishes: a diversity of responses with particular reference to changes in circulating corticosteroids[END_REF]Shreck and Tort 2016). At the physiological level, the primary response to a stressor in fish is neuroendocrinal and involves the activation by the central nervous system of the sympathetic nervous system, releasing neurotransmitters and hormones such as catecholamine and corticosteroids (mainly cortisol in fish). Following the stimulus, plasma cortisol increases after a period of time ranging from a few minutes to an hour, depending on the species and the intensity and/or duration of the stimuli [START_REF] Barton | Stress in fishes: a diversity of responses with particular reference to changes in circulating corticosteroids[END_REF]. Cortisol is considered a good quantitative indicator of acute and chronic physiological stress [START_REF] Ellis | Cortisol and finfish welfare[END_REF]. Until recently, this hormone was mainly measured in the plasma of unstressed and stressed fish, providing a basal and high cortisol level. Yet plasma measures of cortisol can only provide insight into the stress level of the fish at a given time. In an ecological context, an integrative value of stress over a long period is much more valuable, allowing a more comprehensive idea of the life history of the animal. Recently, new matrices have emerged to measure cortisol; these are being used as relatively non-invasive markers of stress [START_REF] Sadoul | Measuring cortisol, the major stress hormone in fishes[END_REF].

Among these markers, fish scales are the most promising matrix as they provide a record of the integrated stress response over time (i.e. from weeks to months, [START_REF] Aerts | Scales tell a story on the stress history of fish[END_REF] and a grid to highlight behavioral ecology regarding exposure to stress factors such as temperature [START_REF] Goikoetxea | Genetic pathways underpinning hormonal stress responses in fish exposed to short-and long-term warm ocean temperatures[END_REF], since the metabolic clearance rate in scales is much slower than in plasma [START_REF] Laberge | Temporal profiles of cortisol accumulation and clearance support scale cortisol content as an indicator of chronic stress in fish[END_REF]. Previous studies have shown that intra-specific differences in basal cortisol are genetically driven (Øverli et al. 2002a) and are generally aligned with differences in behavior or coping style, in which reactive individuals display higher basal cortisol levels compared to proactive individuals [START_REF] Koolhaas | Coping styles in animals: Current status in behavior and stressphysiology[END_REF].

In trout, for instance, long-term cortisol administration was shown to decrease both locomotor activity and aggression (Øverli et al. 2002b). This could be transcribed at the scale level with cortisol accumulation.

Tuna are present in all oceans, and as top predators they play a major ecological role in offshore pelagic ecosystems. Worldwide, tuna species are exploited by both artisanal and industrial fisheries, using different types of fishing gear. Two tuna species, skipjack (SKJ, Katsuwonus pelamis) and yellowfin (YFT, Thunnus albacares), are among the most landed species in the world (FAO 2018). Although the biology and ecology of tuna have been extensively investigated, particularly to provide scientific advice for fishery management and the conservation of stocks, studies measuring cortisol in tuna are scarce and limited to only one species (Atlantic bluefin tuna, Thunnus thynnus) in one environment [START_REF] Kirchhoff | Effects of the first two months of ranching on the health of southern bluefin tuna Thunnus maccoyii[END_REF][START_REF] Addis | Physiologic Responses to Stress and Changes in Atlantic Bluefin Tuna (T. thynnus) Meat Color During Trap Fisheries Capture and Processing in Sardinia (W. Mediterranean)[END_REF][START_REF] Tracey | Capture-induced physiological stress and post-release survival of recreationally caught Southern Bluefin Tuna[END_REF]. By measuring the basal levels of plasma cortisol, these studies showed a rapid accumulation of cortisol in response to an intense one-off event (capture) and focused on the short-term stress reaction to capture. In general, captured fish are exposed to stressful stimuli, e.g. physical trauma and physiological stress [START_REF] Skomal | Evaluating the physiological and physical consequences of capture on post-release survivorship in large pelagic fishes[END_REF], which is worth investigating mainly in an ethical context. Our study chose another line of approach, using, for the first time, fish scales as a matrix to measure cortisol in two wild marine fish, skipjack and yellowfin tuna (male and female), in order to describe the causal effects of size and location on this indicator.

Materials and methods

The sample for the study was based on individuals from two tuna species, YFT and SKJ, captured by handline around anchored fish-aggregating devices (FADs) in the Atlantic (Gulf of Guinea off the Ivory Coast) and the Indian Ocean (off Reunion Island) (Fig. 1A andB) on 10-11 July and 26-29 September 2018, respectively. The total time between the bite on the hook and the killing of a tuna was between 2 and 10 minutes. To minimize the pain, we euthanized the animals as soon as they were boarded, with a blow to the head that immediately killed the tuna (following European directive 2010/63/UE). The number of fish sampled was n= 22 in the Gulf of Guinea (16 YFT, 6 SKJ) and n= 23 off Reunion Island (9 YFT, 14 SKJ).

Onboard, the fork length (FL in cm, from the tip of the upper jaw to the caudal fin fork) and body mass (W in kg) were measured by scientists. Scale samples (under the pectoral fin) were taken by scraping the flank with a blunt scalpel and these were collected in plastic Eppendorf tubes, which were immediately immersed in liquid nitrogen and stored at -80°C until analysis. The handlers wore gloves during the entire process to avoid external contamination.

The cortisol extraction from the scale samples (n=45) was adapted from two different protocols [START_REF] Aerts | Quantification of glucocorticoids in fish scales as biomarkers for chronic stress[END_REF][START_REF] Carbajal | Cortisol detection in fish scales by enzyme immunoassay: Biochemical and methodological validation[END_REF]. Scales were washed by adding 1.5 mL of isopropanol in each tube, agitated with a ceramic ball mill at 10 Hz for 3 min in a Retsch MM 400 tissuelyzer, and centrifuged at 14,000 g and 4°C for 1 min. The supernatant was then discarded. This washing process was repeated three times to fully remove the covering mucus. The samples were then stored at -20°C for 1 h and lyophilized (2 h freezing phase down to -40°C, followed by a rapid increase up to 40°C lyophilization phase in a 1000 mbar vacuum overnight). Once dried, subsamples were aliquoted in a new 2 mL Eppendorf tube and precisely weighed (mean weight and SE: 95.5 ± 9 mg). One ceramic mill ball was then added to each tube and samples were minced to powder in the tissuelyzer for 1 min at 30

Hz. Cortisol extraction was performed by adding 1 mL of methanol to each tube, vortexed for 10 sec and shaken overnight (18 h) at room temperature. Then the samples were centrifuged at 12,000 g for 3 min, and the supernatant transferred to a borosilicate glass evaporation tube. This last step was repeated after adding 1 mL of methanol to ensure maximum cortisol extraction. The contents of the glass tube were then evaporated to dryness by nitrogen flow at 40°C using a MULTIVAP® Nitrogen evaporator (Organomation®). All samples were then reconstituted in 150 μL of 0.01% Bovine Serum Albumine (BSA) in Phosphate Buffered Saline solution (PBS -2.7 M KCl; 0.137 M NaCl; 1.8 mM KH2PO4; Na2HPO4), vortexed 3 times for 15 seconds and stored at 4°C until the ELISA assay. The cortisol was quantified using a competitive ELISA measurement described in [START_REF] Faught | Maternal stress-associated cortisol stimulation may protect embryos from cortisol excess in zebrafish[END_REF] and adapted from a protocol validated for zebrafish (Danio rerio) [START_REF] Yeh | An optimized whole-body cortisol quantification method for assessing stress levels in larval zebrafish[END_REF]. This involved coating 96-well plates with a cortisol monoclonal antibody (1.6 μg.mL -1 ; East Coast Bio, ME, USA) in PBS and storing them at 4°C for 16 hours. The cross-reactivity of the antibody is 0.4% for cortisone, 7.7% for 11-deoxycortisol and 21.9% for corticosterone [START_REF] Sadoul | Measuring cortisol, the major stress hormone in fishes[END_REF]. While this is a high level of cross-reactivity with corticosterone, it should be noted that levels of this hormone have been reported to be relatively low in fish (at least 10 times lower than cortisol, as measured in plasma) [START_REF] Kiilerich | Corticosteroids[END_REF], so the specificity of the monoclonal antibody used here and in other fish studies [START_REF] Yeh | An optimized whole-body cortisol quantification method for assessing stress levels in larval zebrafish[END_REF][START_REF] Faught | Maternal stress-associated cortisol stimulation may protect embryos from cortisol excess in zebrafish[END_REF][START_REF] Best | Maternal cortisol stimulates neurogenesis and affects larval behaviour in zebrafish[END_REF]) should only marginally change results. Standards were obtained by serially diluting hydrocortisone (Sigma) in PBS (0-25 ng.mL -1 ), and 50μL of either standards or samples were added to the wells in duplicate. Each well then received cortisol conjugated with horseradish peroxidase (HRP) and filled with a detection solution that turns blue when in contact with HRP. Sulfuric acid (1N) was used to stop the reaction. Absorption was read at 450 nm (EL800 Universal Microplate Reader, BIO-TEK INSTRUMENTS, INC.). Cortisol ELISA concentrations (in ng.mL -1 ) were normalized by the sample weight of each sample (expressed in ng of cortisol per g of scales). The intra-assay CV for the two plates was below 10%. The standard curve was performed with 0.5, 1, 2, 5, and 10 ng.mL -1 , and the cortisol concentration measured in the scales was correlated with expected cortisol values for plate 1 (r 2 =0.95) and plate 2 (r 2 =0.98).

We first compared cortisol concentration in scales between species caught at the same site and tested for a potential fish-size effect. Then we pooled all the data together to test whether the cortisol concentration in scales varied as a function of size, location or species in the same linear mixed model.

We added the ELISA plate as a random factor, as samples were run on separate ELISA plates. The interaction was removed from the model when not significant. To ensure that all fixed effects were taken into account similarly, we used a type II Analysis of Variance with Satterthwaite's method. All statistical analyses were performed using R software 3.2.2 (R Development Core Team 2013) and package "lme4" for linear mixed models [START_REF] Bates | Fitting Linear Mixed-Effects Models Using lme4[END_REF]. Maps were designed using the "ggmap" package [START_REF] Kahle | ggmap: Spatial Visualization with ggplot2[END_REF].

Results and discussion

We did not detect significant differences in cortisol concentration in scales as a function of fish length for either species (YFT: t=-0.24; p=0.8; SKJ: t=-0.14; p=0.9) (Table 1). Within each site, the cortisol concentration in SKJ scales was higher (although not significantly) than that in YFT scales (p= 0.17 in the Gulf of Guinea and p= 0.09 off Reunion Island), so that when pooling all individuals together in the model accounting for all fixed effects (Cortisol ~ Location + Length + Species + (1|PlateScales)), cortisol concentration was significantly higher (t=1.97; p=0.0558) in SKJ scales compared to YFT scales, independently of their size and location (Fig. 2). However, it is worth noting that this pattern was driven by 4 SKJ individuals between 54 and 59 cm in FL (3 individuals from Reunion Island and 1 individual from the Gulf of Guinea). According to the sampling size of this study, neither the size:species interactions nor the size:location interactions were found to be significant (p=0.5 & 0.6 respectively).

Cortisol in scales has been proven to be an effective marker of chronic stress in various fish species [START_REF] Aerts | Scales tell a story on the stress history of fish[END_REF][START_REF] Carbajal | Comparative assessment of cortisol in plasma, skin mucus and scales as a measure of the hypothalamic-pituitary-interrenal axis activity in fish[END_REF][START_REF] Laberge | Temporal profiles of cortisol accumulation and clearance support scale cortisol content as an indicator of chronic stress in fish[END_REF][START_REF] Hanke | Chronic stress under commercial aquaculture conditions: Scale cortisol to identify and quantify potential stressors in milkfish (Chanos chanos) mariculture[END_REF]). Yet the mechanisms of cortisol clearance in scales are not well described, and to our knowledge only one recent study on goldfish (Carassius auratus) has provided insights regarding its rate [START_REF] Laberge | Temporal profiles of cortisol accumulation and clearance support scale cortisol content as an indicator of chronic stress in fish[END_REF]. That study demonstrated that following a cortisol treatment, plasma levels went back to initial values after only 3 days, while it took 27 days to return to the basal level in scales. This implies that cortisol is somehow catabolized at a much slower rate in scales than in plasma, likely because this steroid is embedded in calcified structures that protect it. The cortisol values found for both species in our study were well within the range of those reported for other fish species, e.g. European seabass (Dicentrarchus labrax) (from 1.7 to 3.4 ng g -1 , Aerts and De Saeger 2017) and goldfish (Carassius auratus) (from 0.25 to 2.10 ng.g -1 , [START_REF] Carbajal | Cortisol detection in fish scales by enzyme immunoassay: Biochemical and methodological validation[END_REF]. Although our study relied on a small dataset, the cortisol concentrations in scales were significantly higher in skipjack (mean: 4.75 ± 6.56 ng.g -1 ) than in yellowfin (mean: 1.65 ± 1.85 ng.g -1 ) for individuals of the same size. Interestingly, the cortisol concentration did not increase with fish size (i.e. age), which argues for a slow -but active -renewal of cortisol in this matrix, as demonstrated by [START_REF] Laberge | Temporal profiles of cortisol accumulation and clearance support scale cortisol content as an indicator of chronic stress in fish[END_REF], rather than an accumulation over time.

Adult skipjack and juvenile yellowfin tuna (the size classes sampled in our study) are known to share many anatomical, physiological and behavioral similarities. They have similar swimming modes, thermal biology, increased rate functions (e.g. standard metabolic rate, aerobic capacity, heart rate and gut clearance) and cardiac physiology [START_REF] Graham | Tuna Comparative physiology[END_REF], which distinguish them from many other teleost species. They have the same vertical habitat utilization [START_REF] Schaefer | Vertical Movements and Habitat Utilization of Skipjack (Katsuwonus pelamis), Yellowfin (Thunnus albacares), and Bigeye (Thunnus obesus) Tunas in the Equatorial Eastern Pacific Ocean, Ascertained Through Archival Tag Data[END_REF]) and usually feed on the same prey species [START_REF] Jaquemet | Do drifting and anchored Fish Aggregating Devices (FADs) similarly influence tuna feeding habits? A case study from the western Indian Ocean[END_REF], although some diet differences have been observed in specific oceanic regions [START_REF] Olson | Chapter Four -Bioenergetics, Trophic Ecology, and Niche Separation of Tunas[END_REF]. They both form large schools and associate with floating objects, but differences in their fine-scale FAD-associative behavior have been observed [START_REF] Rodriguez-Tress | Associative behavior of yellowfin Thunnus albacares, skipjack Katsuwonus pelamis and bigeye tuna Thunnus obesus at fish aggregating devices (FADs) off the coast of Mauritius[END_REF]. Considering the strong similarities between these two species, one might expect similar cortisol metabolism, resulting in similar concentrations of cortisol in scales. Larger datasets are needed to confirm the differences between these two species. In parallel, experiments on captive tuna should be conducted in order to investigate the time evolution of cortisol concentration over several days and weeks. If a species difference in cortisol concentration is confirmed by larger datasets, this may be explained by species-specific cortisol metabolism. This would imply intrinsic differences in the production and/or clearance of cortisol (dependent upon cortisol-binding proteins, target tissue receptors, tissue uptake and catabolism of cortisol) [START_REF] Mommsen | Cortisol in teleosts: dynamics, mechanisms of action, and metabolic regulation[END_REF]). The production/clearance of cortisol and the rate of scale formation (and therefore the associated integration of cortisol) in these two species have not yet been investigated and represent a perspective for future studies.

In our dataset, however, the difference in cortisol concentration observed between the two species was mainly due to four skipjack individuals. Interestingly, comparable high values were also detected in some extremely stressed European seabass (Aerts and de Saeger 2017) and goldfish [START_REF] Laberge | Temporal profiles of cortisol accumulation and clearance support scale cortisol content as an indicator of chronic stress in fish[END_REF]. These high concentrations of cortisol in scales may be the consequence of specific stressful events experienced by these individuals, such as attacks by predators or fasting periods that might have occurred over the preceding days or weeks. Understanding why some individuals of the same species and the same size, in the same fish aggregation around a FAD, have significantly higher concentrations of cortisol (about five times higher) would be valuable to study on a physiological and behavioral level.

Considering that the clearance rate of cortisol in scales is very slow, it may be that individuals with high concentrations of cortisol belong to different schools than the fish with lower cortisol concentrations.

Tuna aggregations at FADs usually consist of several schools [START_REF] Josse | Typology and behaviour of tuna aggregations around fish aggregating devices from acoustic surveys in French Polynesia[END_REF][START_REF] Dagorn | Behaviour of yellowfin (Thunnus albacares) and bigeye (T. obesus) tuna in a network of fish aggregating devices (FADs)[END_REF][START_REF] Moreno | Fish behaviour from fishers' knowledge: The case study of tropical tuna around drifting fish aggregating devices (DFADs)[END_REF]). As fish from the same school are likely to have a relatively long common history (although the stability of a tuna school appears to vary quite a lot, see [START_REF] Bayliff | Integrity of schools of skipjack tuna, Katsuwonus pelamis, in the eastern Pacific Ocean, as determined from tagging data[END_REF][START_REF] Hilborn | Modeling the stability of fish schools: exchange of individual fish between schools of skipjack tuna (Katsuwonus pelamis)[END_REF], it can be assumed that fish from the same school have comparable concentrations of cortisol as they have usually experienced the same events. This hypothesis could be investigated by comparing intraand inter-school variability of cortisol concentrations in individuals of the same species sampled from different free-swimming schools caught by purse seiners (as sampling tuna associated with FADs does not ensure that they are from the same school).

Conclusion

This preliminary study provides the first measurements of cortisol in the scales of wild tuna. While further studies are needed to better understand the mechanisms of scale cortisol clearance and metabolism in the selected species (e.g. a comparison of the cortisol level in scales and in other biological matrices: mucus, plasma, muscle, etc.), this data offers a good start for designing new experiments. In the same way [START_REF] Willis | Application of bioelectrical impedance analysis as a method for estimating composition and metabolic condition of southern bluefin tuna (Thunnus maccoyii) during conventional tagging[END_REF] suggested taking bioelectrical impedance analysis (BIA) measurements to measure the condition of tuna during tagging operations, we suggest sampling a few scales from individuals during conventional or electronic tagging. This would be a quick and harmless operation that should not disrupt tagging operations and would allow measuring the markers of chronic stress (as opposed to the stress due to capture, measured by cortisol in plasma) of tagged individuals. Comparing the large-scale movements (from conventional or archival tags) or the fine-scale movements within arrays of FADs (from acoustic tags) for fish with different concentrations of cortisol would provide key insights into the ecology of these species, particularly the possible links between stress response and behavior. Supplementary Information. Not applicable.

Fig. 1 .

 1 Fig. 1. Tuna species, skipjack (Katsuwonus pelamis) and yellowfin (Thunnus albacares) were sampled in A) the Gulf of Guinea (Atlantic Ocean), 27 nautical miles off Abidjan (Ivory Coast) on July 2018 and B) La Réunion Island (Indian Ocean), 10 nautical miles from shore, on September 2018. Red dots indicate the sampling location of the Fish Aggregating Devices (FAD) where tunas have been captured.

Fig. 2 .

 2 Fig. 2. Cortisol concentrations (ng.g -1 ) measured in the scales of the two tuna species K. pelamis (SKJ) and T. albacares (YFT), in the Atlantic and Indian oceans, correlated to the fish length (in cm).

Table 1 .

 1 Biometrics data (number, length, weight) of the two tuna species (sp), skipjack (SKJ, Katsuwonus pelamis) and yellowfin (YFT, Thunnus albacares)

		sp	nb			ND	length (cm) min-max	length (cm) average ± SD	weight (kg) min-max	weight (kg) average ± SD
	Gulf of Guinea	YFT 16	7	6	3	46.5-66	58.08±5.22	1.6-4.6	3.13±0.87
	Gulf of Guinea	SKJ	6	3	3	0	56-66	59.78±3.95	2.9-4.4	3.52±0.52
	La Réunion Island	YFT 9	3	6	0	61.5-69.6	64.61±2.77	4.1-5.5	4.51±0.45
	La Réunion Island	SKJ 14	4	10	0	41.4-60	51.57±5.62	1.2-4.0	2.45±0.79

Acknowledgments

This study was conducted in the framework of the Tuna MicrobiOME (THE MOME) project, with funding from Montpellier University of Excellence (MUSE) and the Marine Biodiversity and Conservation (MARBEC) research group. The sampling and slaughtering of animals was carried out by professional fishermen working for the French National Research Institute for Sustainable Development (IRD) Exploited Tropical Pelagic Ecosystems Observatory, which is ISO 9001:2015 certified.

The authors would also like to thank the professional fishermen from La Réunion Saint Gilles Fishing Club for their help in fish sampling. We also thank Aurélie Guillou (IRD), who assisted us in sampling some specimens. We are grateful to Antoinette Adingra, Nadège Kouadio and Constance Diaha, colleagues from the Oceanology Research Center of Abidjan, for their contribution in the fish sex assessment.

Finally, many thanks to the editor (Dr. Margaret Docker) and referees, who provided very constructive comments on the initial version of this manuscript, and to Elise Bradbury for the English proofreading.

Funding

The Tuna MicrobiOME project (THE MOME) receives funding from Montpellier University of Excellence (MUSE) and the Marine Biodiversity and Conservation (MARBEC) research group.

Code availability. Not applicable.

Data availability. Not applicable.

Declarations

Ethics approval

All sampling and animal handling respected good scientific practice and animal welfare rules. To minimize suffering, we euthanized the animals as soon as they were boarded. The number of fish sampled was limited to the strict number necessary.

Consent to participate. Not applicable.

Consent for publication. Not applicable.

Conflict of Interest/ competing interests. The authors declare that they have no conflict of interest concerning this manuscript and no competing financial interests that could influence the work reported in this paper.


Fanny Canon

Marie-Bernadette Maillard

Gwénaële Henry

Anne Thierry

Valérie Gagnaire

Positive interactions between lactic acid bacteria promoted by nitrogen-based nutritional dependencies

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

-Positive interactions changed carbohydrate consumption and production of volatile compounds -The study of the nitrogen compounds used by the receiver strains will be further investigated to understand how the proteolytic and non-proteolytic strains positively interact

Figure 3 Figure 5 Figure 1 Figure 2 •

 3512 Figure 3 Three different donor stains resulted in three different types of interactions

  Figure 4


a,b,c,d: statistical differences between cultures in the same row For each property, the color scale from red to green represents the highest and lowest value, respectively.

In a chemically defined medium containing casein and lupin proteins as sole nitrogen sources, positive interactions were observed in co-cultures that associated a proteolytic strain, either Enterococcus faecalis (Efa) or Lactococcus lactis (Lla), and a non-proteolytic strain, Lactiplantibacillus plantarum (Lpl). The positive interactions were mediated by the peptides and amino acids provided by the proteolytic strains and were stronger with Efa (Canon et al. 2021).

The aim of the study was to evaluate the impact of these positive interactions on the organoleptic and textural properties of mixed dairy and plant-based yogurt alternatives (YAs).

As lupin proteins are associated with unpleasant flavour, two different milk:lupin protein ratios were tested (67:33 and 50:50). In order to substitute animal-based fat by plant-based fat, anhydrous milk fat and coconut oil were also tested.

-12 yogurt alternatives were prepared.

-With positive interactions (Efa x Lpl), both strains contributed to the final properties of the YAs. Without positive interactions (Lla x Lpl), only Lla contributed in the final properties of YAs.

-More textured Yas (firmer, less wheying off, more viscous) were obtained with anhydrous milk fat and a milk:lupin protein ratio of 67:33.

-Substitution of milk protein by lupin protein can be acceptable, in terms of flavour and texture, if the milk:lupin protein ratio is about 67:33.

-Coconut oil and anhydrous milk fat gave different physical properties that were not perceived during the sensory analysis. Substitution of animalsourced fat by plant-based one is thus achievable. • Total bacterial counts reached 10 9 cfu/mL in each YA

• Protein ratio and fat type had no impact on bacterial growth (Fig. 1) and acidification (Table 1)

• Lpl reached a significantly higher proportion when co-cultured with Efa compared to Lla (Fig. 1)

• Proteolysis was impacted by the starter and the protein ratio but not by the fat type

• Efa showed a strong proteolytic activity and hydrolysed both milk and lupin proteins (Fig. 2)

• Protein hydrolysis was hardly observed with Lla

• Lpl decreased the levels of free NH 2 compouds

• Only the YAs fermented by Efa did not reach the targeted pH within 12 h incubation.

• Lla acidified faster, followed by Efa and Lpl. In both co-cultures, the acidification rates did not differ from that of Efa and Lla cultures (Table 1).

• Apparent viscosity at 50 s -1 , firmness and WHC were significantly higher for YAs: containing AMF, with a milk:lupin protein ratio of 67:33, and fermented with Lla or Lla x Lpl.

• The association Lpl x Efa resulted in more viscous and firmer YAs, as well as less wheying off, compared to the YAs fermented with Efa alone.

• YAs were differentiated according to the protein ratio (Dim1) and the fat type (Dim2) but not the starter (Fig. 3).

• The protein ratio of 67:33 was described as pleasant, textured (hard gel) and nonhomogenous as opposed to the ratio 50:50 described as unpleasant, bitter and with a mellow texture (Fig. 3).

• AMF was described as milky, lactic and "goaty" and COCO as fruity, fresh and nutty.

Factors Levels Responses

Starter culture

•

  ± 0.01 a 0.29 ± 0.01 ab 0.29 ± 0.01 b 0.27 ± 0.00 ab 0.26 ± 0.01 ab 67:33 0.33 ± 0.01 a 0.44 ± 0.04 b 0.5 ± 0.04 c 0.49 ± 0.05 bc 0.48 ± 0.06 bc

Fig. 3

 3 Fig. 3 Correspondence analysis map of sensory analysis (sorting task) data obtained for the 8 YAs fermented with Lla and Lla x Lpl

Fig 2 .

 2 Fig 2. SDS-PAGE and quantification of free NH 2 groups a,b,c : statistical differences between cultures

Fig 1 .

 1 Fig 1. Percentages of the Lpl strain in co-cultures

  Canon, F., Maillard M-B., Henry G., Thierry, A., & Gagnaire, V. (2021). Positive interactions between lactic acid bacteria promoted by nitrogenbased nutritional dependencies. AEM. doi: 10.1128/AEM.01055-21.

  Enterococcus faecalis CIRM-BIA2412 (Efa) •Lactococcus lactis NCDO2125 (Lla) •Lactiplantibacillus plantarum CIRM-BIA1524 (Lpl) •Co-culture Efa x Lpl •Co-culture Lla x Lpl

Table 1 .

 1 Physical

properties of the 20 yogurt alternatives according to fat, protein ratio and cultures


Jolan Philippe
email: jolan.philippe@imt-atlantique.fr

Massimo Tisi
email: massimo.tisi@imt-atlantique.fr

Hélène Coullon
email: helene.coullon@imt-atlantique.fr

Gerson Sunyé
email: gerson.sunye@ls2n.fr

Executing Certified Model Transformations on Apache Spark

Keywords: CCS Concepts:, Computing methodologies → MapReduce algorithms, • Software and its engineering → Correctness, Domain specific languages Model-Driven Engineering, Model transformation, Spark, Coq

Formal reasoning on model transformation languages allows users to certify model transformations against contracts. Co-qTL includes a specification of a transformation engine in the Coq interactive theorem prover. An executable engine can be automatically extracted from this specification. Transformation contracts are proved by the user against the CoqTL specification and guaranteed to hold on the transformation running on the extracted implementation of CoqTL. The design of the transformation engine specification in CoqTL aims at easing the certification step, but this requirement harms the execution performance of the extracted engine.

In this paper, we aim at providing a scalable distributed implementation of the CoqTL specification. To achieve this objective we proceed in two steps. First, we introduce a refined specification of CoqTL that increases the engine parallelization. We present a mechanized proof of the equivalence with standard CoqTL. Second, we develop a prototype implementation of the refined specification on top of Spark. Finally, by evaluating the performance of a simple case study, we assess the speedup our solution can reach.

Introduction

In model-driven engineering (MDE), model management frameworks propose dedicated languages to transform models, like the AtlanMod Transformation Language [START_REF] Jouault | ATL: A model transformation tool[END_REF] (ATL) or the Epsilon Transformation Language [START_REF] Kolovos | The Epsilon Object Language (EOL)[END_REF] (ETL). Good scalability and facilities for formal reasoning are among the intended benefits for users of model-transformation languages. On the one hand, researchers have proposed transformation engines designed to effectively perform computationally or memory-intensive transformations [START_REF] Dimitris S Kolovos | Scalability: The holy grail of model driven engineering[END_REF]. For instance, transformation languages have been equipped with implicit parallel/distributed modes of execution to automatically multiply the number of resources allocated to a transformation [START_REF] Benelallam | Efficient model partitioning for distributed model transformations[END_REF][START_REF] Kolovos | The Epsilon Object Language (EOL)[END_REF]. To achieve this, some engines have been built on top of distributed programming frameworks (e. g. Apache Hadoop in [START_REF] Benelallam | Efficient model partitioning for distributed model transformations[END_REF]). On the other hand, the community has extensively worked at formal reasoning and verification tools for model transformation languages. Among these solutions, the CoqTL language [START_REF] Tisi | CoqTL: an Internal DSL for Model Transformation in Coq[END_REF] allows users to write transformation rules, define contracts and certify the transformation against them, within the Coq proof assistant [START_REF]The Coq development team[END_REF].

In this paper, we introduce a transformation engine that addresses at the same time distribution and certification. This is not trivial. Distributed programming typically involves non-deterministic execution order, complicating the proof of properties on parallel programs. Interactive theorem proving is already a very costly activity on sequential model transformations. Certifying transformations by formal reasoning on the complexities of modern distributed frameworks would require unmanageable proofs. We propose a solution to allow users to certify transformations in Coq and execute them over a modern data-analytics framework. The core of our proposal is SparkTE, a distributed implementation of the CoqTL specification, built on top of Apache Spark, a widelydeployed framework. Besides scalability, Spark provides an ecosystem of libraries, e. g. for connection to heterogenous data sources. Importantly, the Spark programming interface mirrors a functional/higher-order model of programming. This is a key property for our solution, since it enables a straightforward and high-confidence extraction of parallel code from functional Coq specifications. Notice that we designed the parallelizable specification for a specific data analytic framework, and targeting another back-end solution would necessitate additional manual changes.

This paper has two main contributions:

• A refinement of the standard CoqTL specification, including three major optimizations to increase parallelism opportunities. The refined specification, here named Parallelizable CoqTL, is written in Coq, and is validated by formally proving in Coq the input/output equivalence to the standard CoqTL specification. • A transformation engine that implements Parallelizable CoqTL on top of Apache Spark, named SparkTE. By evaluating the performance of a simple case study, we assess the speedup that SparkTE can reach.

The rest of the paper is organized as follows. We first introduce a running example in Section 2. Our approach overview is illustrated in Section 3. Section 4 presents the new Parallelizable CoqTL specification and the SparkTE engine. The performance experiments on SparkTE are presented in Section 5. We finally conclude the paper by summarizing related work in Section 6 and drawing conclusions in Section 7.

Motivation and background 2.1 Running example

We choose a simple transformation as a running case to illustrate the approach and perform preliminary performance assessments. Listing 1 shows an excerpt of the code of a CoqTL transformation named Relational2Class, which ideally reverse-engineers class diagrams from given relational schemas, and Figure 1 shows its source and the target metamodels. The transformation is a simplified inverse of the well-known Class2Relational transformation [START_REF] Tisi | Parallel Execution of ATL Transformation Rules[END_REF], often used in the community for exemplifying new contributions. We choose the inverse direction because it is representative of reverse-engineering transformations, where scalability problems frequently arise [START_REF] Benelallam | Distributing relational model transformation on MapReduce[END_REF].

The transformation is written in CoqTL, an internal DSL for model transformation within the Coq theorem prover. The transformation primitives are newly-defined keywords (by the notation definition mechanism of Coq), while all expressions are written in Gallina, the functional language used in Coq. The CoqTL semantics is heavily influenced by ATL [START_REF] Jouault | ATL: A model transformation tool[END_REF] (notably in the distinction between a match/instantiate and an apply function), and its original design choices focus on simplifying proof development.

In Listing 1, the Relational2Class transformation is defined via four rules, composed of two parts: (i) a matching section (type and guard condition) and (ii) an output section, which contains a definition for created target elements and optional references. To keep a trace of which expression is used for mapping a source to a target element, each element of the output section of the rules is named.

The first rule (lines 5-11) maps all relational Types to class-diagram Datatypes. We specify that a datatype is constructed using BuildDataType with the same id and the same name as the matched type (line 10). The Table2Class rule (lines 5-11) translates all tables to a corresponding class, with the exception of tables that persist multivalued attributes. To filter these tables, a guard condition, introduced by the when keyword, calls a user-defined isClassTable function (line 14), that we will discuss later. The created target class is constructed from the id and the name of a matched table (line 18). The third rule (lines 21-40) generates Attributes and links them to Classifiers. Columns are transformed into single-valued Attributes by a call to BuildAttribute with the last parameter (multivalued) equal to false (line 27). Again, columns that are contained by tables that persist multi-valued attributes are not transformed, thanks to a guard isClassTable (line 14). Additionally, two links are defined for the generated Attribute. First, Listing 2. The tr_tracePattern_source lemma table t generated from a multi-valued attribute and the table o that corresponds to its owner. The two references are built as in the Column2Attribute rule: one reference to the type of the attribute and one to its owner.

Once the transformation has been defined in CoqTL, it is possible to prove that it has a given property, e. g. it respects a given contract. For instance, CoqTL users can interactively prove a properties of Relational2Class and Class2Relational, by proving the lemma in Listing 2 using Coq tactics. Several strategies are possible to distinguish tables that were generated from classes, from tables that correspond to multivalued attributes. These strategies would correspond to different implementations of the functions isClassTable and isMVAttributeTable. In our simple example, we assume that Class2Relational translates multivalued attributes into tables with a specific name pattern (e.g., name of the owner class followed by an underscore and the name of the multivalued attribute). Under this assumption isClassTable and isMVAttributeTable just contain a string check, with constant complexity. In the experimentation (Section 5) we will introduce more computational time in these functions to simulate increasing complexity, to better estimate the potential parallelization of our solution.

Objective

CoqTL includes an executable semantics for the transformation engine. An implementation of the engine in OCaml or Haskell can be automatically obtained by the standard extraction mechanism of Coq. However, since the executable semantics is designed to be reasoned about in proof terms, it is kept very simple (the core semantics is formalized in 196 Coq LOC) and does not include any efficiency optimization. This has a significant impact on the performance of the extracted engine for CoqTL. For instance , the execution of the Relational2Class transformation on a sequential version of the CoqTL engine for a model of a thousand of model elements takes more than 4 hours on a recent laptop (Intel Core i7-8650U CPU @ 1.90 GHz, 32 GB of RAM).

Considering the size of realistic code-bases for reverseengineering projects, we aim for a solution that allows users to deploy and run the transformation on a state-of-the-art distributed data-analytics framework. Such frameworks are often already deployed in companies that perform large computations, or in their cloud-based services portfolio.

Apache Spark is a data-analytics framework aiming at querying or manipulating large-scale data. In a local mode, the computation is run using a single machine, while in a cluster mode a master dispatches some partitions and tasks to several other additional machines: the workers. Such architecture allows to use more resources (memory, number of processors) which lead to the opportunity of running computations on larger datasets. The data-parallel approach promoted by Spark consists in using a specific data structure called RDDs (Resilient Distributed Dataset) that is partitioned by the master. Partitions are distributed among the workers as well as tasks to execute on RDDs. A task in Spark is a set of operations performed on a RDD. Spark supports several programming paradigms. The use of these paradigm for model-management operations has been investigated in [START_REF] Philippe | Towards Transparent Combination of Model Management Execution Strategies for Low-Code Development Platforms[END_REF]. In this paper, we use Spark with the MapReduce paradigm [START_REF] Lammel | Google's MapReduce programming model -Revisited[END_REF].

We aim at designing a solution that executes on Spark a CoqTL transformation (like Relational2Class), certified against a contract. On the performance side, we aim at improving two kinds of scalability: 1) The capacity of scaling with additional computational resources, referred to as vertical scalability; 2) The capacity of dealing with increasing datasets, referred to as horizontal scalability. On the reliability side, while the whole execution stack can not be certified end-to-end (e. g. the formal semantics of Apache Spark is not available), we aim nonetheless to produce a solution that gives users high confidence that the proved contract will hold on the distributed transformation.

Approach overview

Figure 2 shows the global workflow of our approach, separating the formal part, written in Coq (left side), from the executable engine in Scala and Spark (right side). Starting from the Coq part, the workflow must be interpreted as follows.

The CoqTL language allows users to define model transformations, theorems on their behavior and machine-checked proofs of these theorems in Coq. These transformations can be executed directly on the executable CoqTL specification. We introduce Parallelizable CoqTL, a refinement of the executable CoqTL specification that increases its parallelization opportunities. We use the Coq theorem prover to prove that this new specification is a refinement of the original one (see Section 4.2). This entails that, for a given source model, any transformation produces the same output model when it runs on standard CoqTL or Parallelizable CoqTL. The Parallelizable CoqTL specification is written in Gallina, the functional language used in Coq.

Then we provide an implementation of Parallelizable Co-qTL in Scala. To obtain the Scala implementation, that we name ScalaTE, we manually extract the Gallina functions into corresponding Scala functions. Section 3.1 details this extraction.

ScalaTE uses data structures that are the direct correspondent of Gallina data structures (e. g. Scala List for Gallina list). In a final step we replace these data structures with distributed data structures (RDDs) from the Spark library, as described in Section 3.2. We name the resulting distributed engine SparkTE.

To run a transformation on SparkTE, the CoqTL transformation rules (e. g., Listing 1) have to be translated into their corresponding Scala version. We currently perform this step manually, but its automatization is possible and left for future work. Since Spark does not change the functional interface of standard Scala, the obtained Scala transformation can run on both ScalaTE and SparkTE.

Coq to Scala

The Coq environment includes an extraction mechanism targeting ML languages: OCaml, Haskell, or Scheme. Although an automatic extractor to Scala is available [START_REF] El | Scallina: Translating Verified Programs from Coq to Scala[END_REF], it supports only a subset of Gallina on an outdated version of Coq. Hence, we opted to perform the extraction manually. We perform manual extraction at two levels: first to create the core engine (gallina2scala in Fig. 2), then to obtain Scala rules representing a CoqTL transformation (coqtl2scala).

Gallina to Scala. The executable CoqTL specification can be seen as a functional program that interprets the transformation code. We produce a literal translation of this interpreter in Scala.

For extracting Scala code from Parallelizable CoqTL we translate Gallina types and (pure) functions into their correspondent types and pure functions in Scala. Listings 3 and 4 show an example of CoqTL executable specification (in Gallina) and its implementation in SparkTE (in Scala). The example well illustrates the one-to-one translation among types and functions that we adopt for all the extraction.

We implicitly assume that Scala functions (e. g., flatMap for Scala List) implement the same semantics than their Gallina correspondent (e. g., flat_map for Gallina list).

CoqTL to Scala. The CoqTL parser translates the concrete syntax of the transformation (e. g., from Listing 1) into Coq code to construct an abstract syntax tree. Obtaining the same transformation in Scala requires constructing the same abstract syntax tree as Scala objects. Note that Scala constructors for the abstract syntax are the literal translation of the corresponding Gallina constructors in CoqTL.

Listing 5 shows the translation of the Table2Class rule from Listing 1. The rule constructor (RuleImpl) requires a name, a list of types for an input pattern (the types argument), a guard condition (the from function), a list of output pattern elements (the to argument). Each output pattern element has a name, and a function for creating output elements from input pattern elements. While not shown in this example, it can be also accompanied by a list of functions for creating links in the output model. As illustrated by the example, the only non-trivial part of this extraction is the translation of the body of expressions for guards and output element creation (e. g. the body of the

Listing 3. CoqTL definition for generating all combination of elements

def tuples_of_length_n[ A] ( s1: List[ A], n: Int) : List[ List[ A]] = n match { case 0 => List( List()) case n1 =>
prod_cons( s1, tuples_of_length_n( s1, n1-1)) } Listing 4. Extracted Scala code from the Coq function tuples_of_length_n anonymous functions in Listing 5). An automatic compiler from CoqTL rules to Scala is a part of future work (Section 7).

Distributed data structures

Spark RDDs are data structures that are automatically partitionned and resulting in the distribution of the computation operations on a Scala sequence (e. g., List) of serializable elements. From a user point of view, RDDs can be manipulated as lists, using the same primitive functions, and parallelism is implicit. The advantage of using such abstraction for parallelism is the semantics preservation of the operations on the distributed structures. Because of the popularity of Spark and its support, we assume the correctness of parallel operations on the data.

An efficient use of RDDs requires an effective partitioning of data. For instance, to take advantage of the internal multithreading mechanism, it is typically recommended in Spark to assign four data partitions to each core. Each independent computation of a partition is referred as a task. The Spark Listing 5. Scala implementation for the Table2Class rule task scheduler makes the distribution following a roundrobin approach, optimizing load-balancing: once a task is ended on a core, a new one can be assigned from the waiting list. Section 4 gives more detail about how we use RDDs in SparkTE.

Parallelizable semantics for CoqTL

The execute function shown in Listing 6 is the entry point of the transformation execution in the standard CoqTL semantics.

First, the allTuples function (line 3) produces all the tuples of elements that can be possibly matched by the rules. allTuples computes a list of 𝐴 𝑎=0 𝑛 𝑎 tuples, with 𝑛 the number of elements in the input model, and 𝐴 the maximum arity of the transformation rules.

Then, the instantiatePattern function (line 5) tests each tuple to find if it matches with any rule and for each match it constructs the corresponding output pattern elements. Internally it iterates on each rule, executes their guard function and if the result is positive, executes the element creation function for each output pattern element of that Finally, the applyPattern (line 6) function is executed on each tuple to create target links. Similarly to the function instantiatePattern, the function internally iterates on all rules and checks if the rule matches the given pattern. In a positive case, the element creation functions for that rule are executed and then the link creation functions. The resulting links are gathered by the flat_map in a single output list.

Parallelizable CoqTL Parallelizable CoqTL contains three optimizations to the base CoqTL specification:

• To increase parallelization, the algorithm is split into two consecutive phases, instantiate and apply, that are built on parallelizable functional patterns (flat_map);

• To improve load balancing of the instantiate phase, only possibly useful tuples are generated and then distributed; • To improve load balancing of the apply phase, a set of trace links is produced by the instantiate phase and the apply phase iterates only on those trace links.

Note that similar optimizations (among others) are already implemented in well-known transformation engines, like ATL [START_REF] Jouault | ATL: A model transformation tool[END_REF] or ETL [START_REF] Kolovos | The Epsilon Object Language (EOL)[END_REF]. Differently from previous work, in this paper we formalize the optimizations, interactively prove that they do not affect the transformation output and assess their impact on distributed execution.

The entry point of Parallelizable CoqTL, is presented in Listing 7. Figure 3 illustrates the global behavior on a minimal example with one table and one column.

Optimization 1: Two phases. In standard CoqTL the applyPattern function performs all the computation of the links generated by a matched input pattern, by the rule that matches it. However the computation of the links of a rule is not independent from the computation of other rules. This dependency is caused by the resolve function (e. g., lines 32 and 37 in Listing 1) that searches for the output of another rule in order to set the target of the created link. In general, because of this dependency, two executions of the apply function can not be run in parallel, without replicating some matching and instantiation within each call to resolve. 1 We refactor the computation to split it in two phases, similarly to ATL [START_REF] Jouault | ATL: A model transformation tool[END_REF]. This is visible in Listing 7. In the first phase (lines 5) we compute the tuples and we run the matching and instantiation by a new function named tracePattern. The first phase produces the list of generated elements, and trace links connecting them to their corresponding source patterns. Differently from Listing 6, here the second phase 1 This is exactly how the standard CoqTL specification computes resolve. This simple solution keeps the specification compact and has no negative impact on proofs. However it has a big impact on performance and parallelization.

(line 8)) can only start computing output links after the full first phase has finished computing the trace links, since the flat_map expects the tls structure as parameter.

In Listing 7 every execution of tracePattern can be run in parallel. When the first phase is over, every execution of applyPatternTraces can be run in parallel too, since the calls to resolve can be computed immediately on the tracelink structure. This greatly improves the parallelization of the algorithm.

Optimization 2: Tuple generation by rule. Matching a pattern to a rule happens in two consecutive steps. First, the types of the pattern are checked against the types expected by the rule. Then, if the types are correct, a guard condition is evaluated. The type checking is very fast, hence it acts as a first filtering. Instead the evaluation of the guard condition can potentially be very long, or navigate large parts of the model. So it is executed only for the few tuples that pass the type check.

Since the tuples that require an evaluation of the guard condition are a small subset of all the possible tuples, arbitrarily distributing all tuples among the cores can potentially lead to imbalanced partitions. In particular it would not be uncommon to have partitions that do not require any guard evaluation, opposed to partitions that need to evaluate several expensive guards. In such cases, idle workers would wait for the synchronization barrier of Spark to start new computations. The imbalance impacts the scalability of the program.

To limit imbalance, in the initial sequential tuple generation phase, we generate only tuples whose type matches with at least one rule of the transformation. This is shown in Listing 7 by the use of the allTuplesByRule for tuple generation (line 2). allTuplesByRule iterates on rules and produces only combinations of elements of the types listed in the rule input pattern. This improves load balancing of the first phase since all the produced tuples require a guard evaluation.

Optimization 3: Apply iterates on traces. Executing the apply phase on the tuples generated by allTuplesByRule would cause an imbalance among partitions, similar to the one discussed in Optimization 2. Indeed, among these tuples only very few have passed the guard condition in the first step. A partitioning of allTuplesByRule would produce partitions that do not require any computation, together with partitions that need to evaluate several expensive link creation functions. In this optimization we make the apply phase iterate only over the source patterns that passed the guard evaluation in the instantiate phase. We retrieve these patterns by collecting them from the list of trace links. This is performed by the function allSourcePatterns at line 8 of Listing 7.

Refinement proof

Besides the executable functional specification, CoqTL is also described by an axiomatic specification. Certifying against the axiomatic specification involves providing 10 types, 27 semantics functions and proving 15 theorems. The specification is fully illustrated in [START_REF] Cheng | Certifying a Rule-Based Model Transformation Engine for Proof Preservation[END_REF], together with a proof that engines implementing the executable specification certify against the axiomatic one.

We prove that engines implementing Parallelizable Co-qTL certify the axiomatic specification, too. For this step, we naturally reuse the types, functions and certification proofs of the base executable specification that are not impacted by the optimization. Each optimization is proved independently. Table 1 shows the size (in lines of code) of new specifications and proofs required for describing and certifying each optimization, plus the human effort (in man-days) to complete the proofs. The refined specifications and their proofs are Listing 8. Certifying optimizations lemmas available online 2 . All the discussed lemmas below can be found in Listing 8.

The key step for certifying the twoPhases optimization, is proving that it does not change the output of the full transformation (exe_preserv). This proof has two parts: 1) for the instantiation phase, we prove that the additional computation of the trace, does not have any effect on the computation of the instantiated elements; 2) for the apply phase we need to prove that the new apply function is equal to the old one. Coq is able to prove automatically the second step, by full unfolding and simplification of the old and new apply functions. Note that in this proof we use the axiom of functional extensionality (two functions are equal if their values are equal at every argument), notably to compare the body of inner anonymous functions.

The byRule optimization changes the order of the element and link creation in the transformation output, hence a lemma similar to exe_preserv would not hold. We prove its equivalence in two steps: 1) we prove that tuples computed by rule are a subset of the all tuples (In_by_rule), 2) we prove that tuples computed by rule include all tuples that produce elements or links (In_by_rule_instantiate). The second step is the most challenging and is performed by case analysis on matchPattern: given any pattern, if the pattern does not match then it cannot produce anything, if it matches then it is one of the patterns generated by the rule allTuplesByRule.

The iterateTraces optimization does not change any function in the instantiate phase, hence we can easily prove a similar theorem to (exe_preserv), but limited to the result of instantiate. The apply phase does not change the computation of a single link, but skips some source patterns, and produces links in a different order w.r.t. the base version. To prove that the same links are generated, we proceed by case analysis on the trace generated for a given source pattern: 1) if the instantiate phase did not generate any trace for that pattern, then the standard apply phase would not have produced any corresponding link; 2) if there is a trace, then the apply phase applies the same function to the same source pattern, hence producing the same output (tr_tracePattern_source).

Implementation on Spark

Distributed computation of our Spark implementation revolves around the use of RDDs. First, the input patterns, represented as a list of tuples, are distributed with an RDD among cores. Each core independently applies the instantiate function to every tuple of its partition. Implicit communications are operated by Spark to scatter the tuples from the master process to the workers. After the computation, the resulting elements and their trace-links are all gathered to the master process. Table 3. Description of the three datasets used in the experiments with the number of elements and links.

For the second phase, the trace-links are distributed with an RDD. Each core is in charge of generating the links for a partition of output elements and their associated applied rule. Since this second phase needs a global knowledge of the trace-links to resolve output elements, we used a broadcast communication to share the whole set of trace-links to all cores.

A global view of where RDDs are used and what communications are operated is illustrated in Figure 3.

Limitations

SparkTE implements the complete CoqTL specification but some limitations remain. First, all the translations from Coq to Scala are manual. In particular, the transformation rules need to be manually translated by the user. However the correspondence between the abstract syntax in CoqTL and ScalaTE/SparkTE makes the translation of the rule structure trivial. Hence, our solution reduces the complex problem of translating to Scala a CoqTL transformation, into a simpler problem: translating to Scala separately each side-effect-free expression for guards, output pattern elements creation, and output pattern links creation. While simpler, this task is still tedious and error prone, hence we aim at making the translation fully automatic in future work.

Experiments

In the previous sections, we have presented the three optimizations adopted for Parallelizable CoqTL to increase the efficiency and propensity to parallelization of the extracted engine. In this section, we evaluate the performance of SparkTE and the quality of its speedup compared to the ideal speedup that divides the execution time by the number of cores (relative to a reference version). As none of the contributions of the related work have studied the use of Spark for parallel transformations, we restrict our comparisons to this theoretical ideal speedup. All jobs are run with Spark in a standalone mode.

The first subsection presents the results we get from the use-case previously described in Section 2.1, i. e. the Rela-tional2Class transformation , and similar experiments on two additional cases.

We show that the three optimizations introduced in Section 4 improve significantly the performance of the transformation in Spark. The results also show a relatively poor speedup compared to the ideal one, but that can be improved when having more complex operations within the transformation. To investigate this result, the second subsection illustrates that a speedup close to the ideal can be observed with a high number of cores if the computation time and the size of the dataset is big enough to offset the overheads of Spark. Finally, a performance analysis by phase is presented to show the very good speedup obtained by the parallelized steps of SparkTE.

All results presented in this section have been executed ten times on the same hardware configuration, and an average is presented. The standard deviation of our measurements are never higher than 10 % of the total execution time, thus guaranteeing our results to be stable, and the average to be meaningful. Furthermore, all our experiments have been conducted on the French experimental platform for distributed computing Grid'5000. Grid'5000 is made of geographically distributed clusters of machines in France, each one with its specific hardware setup. This platform also facilitates the reproducibility of the experiments by offering a way to build the same environment for any researcher. However, Grid'5000 makes us dependent on available machines (i. e. nodes) for provisioning which is why we use two different clusters of Grid'5000 in our experiments. For each experiment the number of machines and the number of cores are specified. One can note that the number of machines that is specified correspond to the number of workers instantiated in Spark and that one additional dedicated machine is provisioned to host the master of Spark. All our codes, raw results, and analysis scripts are publicly available online. 3

Evaluation of SparkTE on use-cases

In this section we apply our running example to a first model, denoted D1 in Table 3. We also consider two additional transformations: the IMDb case [START_REF] Horn | The TTC 2014 Movie Database Case[END_REF], aiming at finding couples of actors who recurrently played together, and queries on a DBLP model to find active authors who published in specific journals [START_REF]Lintra[END_REF]. We run the two additional transformations on model instances from the IMDb and DBLP metamodel, respectively denoted D4 and D5 in Table 3.

The transformations specified in CoqTL, subsequently translated to Scala and Spark, are always correctly computed by SparkTE, i. e. with the expected output of the transformation, for all our experiments.

We have used from 1 to 2 machines of the cluster paravance detailed in Table 2. Table 4 presents the execution times and speedups of the three transformations from 1 to 8 cores with 1 and 2 machines. The results show a relatively poor speedup. It is caused by the lack of complexity in the operations used in the guard condition, the instantiate function and the apply function.

Spark's overheads are more particularly observable with the poor speedup measured with 2 cores. Indeed, in this case the gain of parallelizing with Spark is completely swallowed by the overheads, but then an improvement is obtained with 4 cores. Overall the obtained speedup is not satisfying as the speedup for 8 cores is almost the same as with 4 cores, and far from the expected ideal speedup of 8. The additional cases show that the performance of the approach depends on the computation time of the rules. In the IMDb case, the better growing speedup is caused both by the rules themselves that have a high-level of complexity, and the high connectivity of the input model. On the contrary, in the case of the DBLP transformation, the rules with a low computational cost lead to a poor speed-up. We will show in the next subsection that Spark's overhead is fully compensated by high computation time of the transformation and/or large size of the dataset.

One problem with computation on large models is the lack of memory. Distributed solutions, such as the one we propose here, allow users to increase the resources allocated to a transformation. For instance, a fixed allocated memory for one core is not enough for processing the DBLP example, causing disk operations which drastically slow the computation. This phenomenon disappears with 2 cores, and we observe a hyper speedup.

Finally, we want to show the performance gain obtained thanks to our three optimizations. To this purpose, we execute the same naive transformation on the direct implementation of the CoqTL specification (without optimizations). On 1 core naive is computed in 27 s by SparkTE and in 52 s by the CoqTL implementation.

Performance analysis by complexity and datasets

In the previous subsection, we have shown that the performance of SparkTE is enhanced by our three optimizations in terms of execution time and speedup. However, compared to the ideal speedup our use-cases are disappointing. In this subsection, we demonstrate that these results are due to the small computation time in the transformation and the small size of the dataset.

Experimental setup. To this purpose, we have built a version of the transformation with additional fictitious processing time in the different phases and steps. More precisely, we incorporate sleeping functions to different parts of the q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 10 20 30 1 2 4 8 16 32

number of cores relative speedup

Introduced simulated computation time on phases q q q q q q q 0ms 50ms 120ms 250ms 500ms 1000ms 2000ms ideal speedup 50 of ideal speedup (a) Relative speedups of SparkTE for D1 up to 32 cores on 4 machines q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 10 20 30 4 8 16 32 64 128 number of cores relative speedup (b) Relative speedups of SparkTE for D2 up to 128 cores on 8 machines q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 10 Experiments respectively conducted with the dataset D1, D4, and D5, on the cluster paravance.

rule evaluation: the guard condition to simulate complex checking functions; the instantiate part to simulate complex instantiation of element; and the apply function to simulate a long resolution for links. As we discussed in Section 2, userdefined functions can have several implementations with different complexities. Proposing an evaluation based on the computation time of these functions offers an infinite set of benchmarks to accurately estimate the speedup of SparkTE. Furthermore, we additionally vary the size of the dataset in the obtained benchmarks. Table 5 summarizes the set of benchmarks explained hereafter. The first benchmark (B1) uses the same data introduced in the previous subsection on (D1) with 1, 2, and 4 machines, and up to 8 cores per machine on the paravance cluster of Grid'5000. We compare the computation time of the transformation on 2, 4, 8, 16, and 32 cores relative to the execution time on 1 core. The second and third benchmarks, (B2) and (B3), illustrate both the horizontal scalability and the vertical scalability of SparkTE by increasing the size of the dataset with (D2) and (D3) detailed in Table 3, and the number of machines (i. e. nodes, workers) and cores detailed in Table 5. The two benchmarks are evaluated on 1 to 4 machines, up to 16 cores per machine. As increasing the size of the dataset and the fictitious sleeping time also increases the total execution time of experiments, we compute our speedup relative to 4 cores instead of 1. As a result, the ideal speedup for 16 cores, for instance, is 4. Experiments on (B2) and (B3) are conducted on the gros cluster of Grid'5000. Furthermore, as indicated in Table 5, the three benchmarks are evaluated on the following sleeping times: 0, 50, 100, 250, 500, 1000, and 2000 ms. Finally, for the three benchmarks the overall execution time is measured, thus including some small sequential parts of the code, and the broadcast phase.

Results. Figures 4a,4b, and 4c show the speedups observed for each benchmark according to the sleeping time and compared to the theoretical ideal speedup. When increasing the sleeping time, i. e. the execution time of the transformation, the speedup is enhanced and is closer to 5. The set of experiments described in Section 5.2 with different fictitious processing time in the transformation, and different sizes of datasets. The number of cores used by each benchmark is also indicated, as well as the number of cores per node (i. e. worker). Finally, the Grid'5000 cluster used for the benchmark is given in the last column.

the ideal one. By increasing the size of the dataset, one can note a slight increase of the speedup by comparing B2 and B3 results respectively in Figure 4b and Figure 4c. Indeed, in Figure 4b, at 128 cores, more than half of the points are below the 50 % optimal value, while in Figure 4c only two points are below this theoretical value. In other words, at 128 cores, (D2) needs a sleeping time of 500 ms to reach 50 % from optimal speedup while (D3) only needs a sleeping time of 120 ms.

Performance analysis by phase

This third subsection aims at analyzing the impact of the two phases on the global speedup of the program. To do so, we processed the same benchmark as before (i. e., (B1), (B2) and (B3)), but with additional counters within the program to record the computation time of each distinct part: (1) the tuples generation; (2) the instantiate phase; (3) the broadcast intermediate step; and (4) the apply phase.

Results. Our results on the three benchmarks show that the total computation time is mainly driven by the parallel parts. Indeed, in the case of the biggest dataset, with a sleeping time of 0 ms, the sequential part represents only 3.5 % of the total computation time which is the maximum percentage of the sequential parts. For instance, at the opposite extreme, i. e. in the case of the smallest dataset, with a sleeping time of 2000 ms, the sequential part represents less than 1 % of the total execution time. Hence, in the following, we restrict our analysis to the speedup of the parallel parts.

Table 6 illustrates the impact of each phase, by comparing their relative speedup to the optimal one. For reading convenience, we only show the results for sleeping times equal to 50 ms and 2000 ms. Table 6 confirms our previous assumption about the overhead of Spark that can be offset by increasing the computation time (i. e. sleeping duration) or the size of the dataset. Also, one can note that the apply phase offers better scalability than the instantiate phase. This result can be explained by the remaining imbalance on the tuples distribution. Let us remind first that the first instantiate phase is composed of two steps: (1) a guard condition; (2) the instantiation of the tuples that have satisfied this condition. Hence, even if the partitions are better balanced by our second optimization (see Section 4.1), not all the tuples are necessarily computed. As a matter of fact, the guard condition is always evaluated, but not the instantiation that depends on the result of the guard condition. On the contrary, in the apply phase apply patterns to all entries which are nearly perfectly balanced.

Related work

Since models are basically typed graphs, several research efforts rely on graph theory, especially for transformations. For instance, in their Visual Modeling and Model Transformation (VMTS) tool, Mezei et al. use graph rewriting operations based on task-parallelism to apply distribute matching operations to large models [START_REF] Mezei | Towards truly parallel model transformations : A distributed pattern matching approach[END_REF]. Recently, the Pregel [START_REF] Malewicz | Pregel: A System for Large-scale Graph Processing[END_REF] paradigm, a strategy that aims at easing parallel computations on graphs on top of MapReduce architectures, has been used as a solution for model transformation. This is the case of Krause et al., who describe how transformation rules from the Henshin framework can be partially extracted to Pregel code [START_REF] Krause | Implementing Graph Transformations in the Bulk Synchronous Parallel Model[END_REF]. A second use of Pregel in model transformation is proposed by Thung and Hu [START_REF] Tung | Towards Systematic Parallelization of Graph Transformations Over Pregel[END_REF]. In their approach, graph transformations specified in a DSL are transformed into executable Pregel code.

Benelallam et al. introduce ATL-MR, a tool that uses the MapReduce paradigm to distribute the execution of ATL transformation rules [START_REF] Benelallam | Distributing relational model transformation on MapReduce[END_REF]. Their data-parallel implementation follows a similar strategy as our SparkTE implementation: sub-parts of the model are independently transformed in a map phase and then dependencies are resolved in a reducing phase. In addition to distributing computation, authors have recently highlighted the good impact of their strategy for data partitioning. For instance, Burgueño et al. propose LinTra, a prototype for transformations on distributed architectures that tackles memory issues they face with sharedmemory solutions [START_REF] Burgueño | Towards Distributed Model Transformations with LinTra[END_REF]. They applied different strategies mixing both the a distribution of tasks and data, on a single or on multiple machines.

Both, ATL-MR and LinTra show better performances than SparkTE for executing a model transformation. However, they do not propose a verification counter part for certifying their distributed transformations.

Interactive mechanized proof for model transformations properties is an active research area. Calegari 6. Relative speedups of SparkTE parallel phases on B1, B2 and B3 for sleeping times equals to 50 ms and 2000 ms. The percentage of the observed speedup compared to the theoretical ideal speedup is indicated for each result (higher the better).

ATL model transformations and OCL contracts into Coq to interactively verify whether the transformation is able to produce target models that satisfy the given OCL contracts [START_REF] Calegari | A Type-Theoretic Framework for Certified Model Transformations[END_REF]. Stenzel et al. propose a Hoare-style calculus, developed in the KIV prover, to analyze transformations expressed in (a subset of) QVT Operational [START_REF] Stenzel | Formal verification of QVT transformations for code generation[END_REF]. UML-RSDS is a tool-set for developing correct-by-construction model transformations [START_REF] Lano | A framework for model transformation verification[END_REF]. It chooses well-accepted concepts in MDE to make their approach more accessible to model transformation developers. Once the development is achieved, transformations are verified against contracts by translating both into interactive theorem provers. None of these research efforts addresses proving the equivalence of the sequential and the distributed executions of a transformation.

Finally, one can argue that we could have chosen any other back-end language or framework instead of Spark. In particular, some back-ends automatically extracted from Coq, which would enhance the automatic certification of our pipeline. For instance, we could have chosen to extract Haskell code from CoqTL and then use the Haskell distributed parallel Haskell (HdpH) language to introduce parallelism and distribution. Similarly, we could have extracted OCaml code from the CoqTL specifications and then use the BSML [START_REF] Loulergue | Bulk Synchronous Parallel ML: Modular Implementation and Performance Prediction[END_REF] library for parallelization. This is indeed possible and could be the subject of future work. Note that the three optimizations introduced in this paper, and proven equivalent to the initial CoqTL specification, would be useful for any back-end that introduces parallelism and distribution. However, we want to highlight that the goal of our work is to bridge the gap between certified model transformations and data analytics frameworks, such as Spark. Indeed, as explained in the introduction, choosing Spark paves the way towards largescale model transformations. Furthermore, Spark is broadly adopted by companies and is available on most Cloud platforms.

Conclusion and future work

In this paper, we presented a refinement of the CoqTL specification, designed for optimizing the parallel execution of model transformations on Spark. We have illustrated the benefits and the scalability of our proposed optimizations through the Relational2Class example.

In future work, we plan to continue experiments with other use cases (e. g., the IMDB case study). We plan to write a compiler from CoqTL to Scala, to automatically obtain an executable transformation from a Coq specification. Finally, we will study other optimizations, leveraging the vertexcentric paradigm supported by Spark (i. e., GraphX), and the integration with persistence solution (e. g., HDFS).

Figure 1 .

 1 Figure 1. Relational and Class diagram metamodels

  c class Column when 23 isClassTable (getOwner c m) 24 to [ 25 output "svattr" 26 element a class Attribute := BuildAttribute 27 (getColumnId c) (getColumnName c) false ;

Figure 2 .

 2 Figure 2. Global overview of our workflow to execute certified model transformations on Apache Spark.

Listing 6 .

 6 execute function in the base CoqTL specification Definition execute ( tr: Transformation) ( sm : SourceModel) ( sm : SourceMetamodel) : TargetModel := let tuples := allTuplesByRule tr sm mm in let ( elements, tls) := flat_map ( tracePattern tr sm mm) tuples in let links := flat_map ( fun sp ⇒ applyPatternTraces tr sm sp tls) ( allSourcesPattern tls) in Build_Model elements links. Listing 7. execute function in the Parallelizable CoqTL specification rule. The resulting elements are gathered by the flat_map in a single output list.

  Lemma exe_preserv : ∀ ( tr: Transformation) ( sm : SourceModel), twophases. execute tr sm = execute tr sm. Lemma In_by_rule : ∀ ( sp: list SourceModelElement) ( tr: Transformation) ( sm: SourceModel), In sp ( allTuplesByRule tr sm) → In sp ( allTuples tr sm). Lemma In_by_rule_instantiate : ∀ ( sp: list SourceModelElement) ( tr: Transformation) ( sm: SourceModel), In sp ( allTuples tr sm) ∧ instantiatePattern tr sm sp <> nil → In sp ( allTuplesByRule tr sm). Lemma In_by_rule_apply : ... Lemma tr_tracePattern_source: ∀ ( tr: Transformation) ( sm : SourceModel) ( tl : TraceLink) ( sp: list SourceModelElement), In tl ( tracePattern tr sm sp) → sp = TraceLink_getSourcePattern tl.

  Table; element o class Table Excerpt of the Relational2Class Transformation in CoqTL a link from the target Attribute to its type (lines 29-33).The corresponding type is defined in BuildAttributeType by resolving the output of the Type of the source Column (line 32). The second reference, from the Attribute to its owner (lines 34-38), is defined in the same manner (i. e., by resolving the output of the owner Table of the source Column) (line 37). MVAttributeTable2Attribute (line 42-49), that is the last rule of our transformation, matches two tables: a

		Lemma tr_r2c_inverse: ∀ ( m : ClassModel),
		execute Relational2Class ( execute Class2Relational m)
		= m.
	45	when isMVAttributeTable t o
	46	to [
	47	output "mvattr"

48 element a class Attribute := ... 49 ] Listing 1.

  Distributed computation of a model transformation on SparkTL. 𝑐 denotes a column, 𝑡 a table, 𝑎 an attribute, and 𝑐𝑙 a class. Green elements and links illustrate the entities created by the transformation in the output.

			RDD		RDD			RDD
	t u p l e s	t t	c	instantiate + trace-links	t	cl	apply	t	cl
	t	t	t					c	a	cl
						broadcast
	c	c						t	cl	a
	t u p l e s	c c	c t	instantiate + trace-links	c	a	apply	c	a
	Figure 3.

Definition execute ( tr: Transformation) ( sm: SourceModel) : TargetModel := let tuples := allTuples tr sm in let elements := flat_map ( instantiatePattern tr sm) tuples in let links := flat_map ( applyPattern tr sm) tuples in Build_Model elements links.

Table 1 .

 1 Size

	Optimization Spec. size Cert. size Proof effort
	twoPhases	69	484	10
	byRule	42	487	7
	iterateTraces	69	520	4
	(in LOC) of new specification and certification
	proofs added for each optimization, with proof effort (in
	man-days).

Table 2 .

 2 Hardware setup of the two clusters used during experiments. These clusters are part of the Grid'5000 experimental platform for distributed computing.

	Cluster		CPU	RAM	Netw.
	Rennes		2× 8 cores/CPU	128GB	2×
	paravance Intel Xeon E5-2620		10Gbps
	Nancy		2× 18 cores/CPU	96 GB	2×
	gros	Intel Xeon Gold 5220		25Gbps
	Dataset		D1 D2 D3	D4	D5
	model type	Relational	IMDb DBLP
	elements		150 300 600	440	700
	links		290 580 1060 1968	1886

2 https://github.com/atlanmod/SparkTE_public/

Table 4 .

 4 Execution times and speedups (relative to 1 core) for the executions with SparkTE of Relational2Class (R2C), the IMDb findcouple transformation, and the DBLP case.

		30
	relative speedup	20
		4 8 16	32	64	128
				number of cores
	(c) Relative speedups of SparkTE for D3 up to 128 cores on 8 machines
	Figure 4. Relative speedups of SparkTE with sleeping times

Table

		dataset	sleeping	♯ cores	♯ machines ♯ cores per node	cluster
	B1	D1	{0, 50, 120, 250, 500, 1000, 2000}	{1, 2, 4, 8, 16, 32}	{1, 2, 4}	4 to 8	paravance
	B2	D2	same as B1	{4, 8, 16, 32, 64, 128} {1, 2, 4, 8}	4 to 16	gros
	B3	D3	same as B1	same as B2	same as B2	same as B2	gros

  et al. encode Cores Instantiate phase (50ms) Apply phase (50ms) Instantiate phase (2000ms) Apply phase (2000ms)

		1	1 (100%)	1 (100%)	1 (100%)	1 (100%)
	B1	2 4	1.48 (74%) 2.40 (60%)	1.62 (81%) 1.83 (45.75%)	1.81 (90.5%) 3.21 (80.25%)	1.87 (93.5%) 3.85 (96.25%)
		8	3.40 (42.5%)	4.50 (56.25%)	5.78 (72.25%)	7.15 (89.375%)
		4	1 (100%)	1 (100%)	1 (100%)	1 (100%)
		8	1.68 (84%)	1.90 (95%)	1.83 (91.5%)	1.99 (99.5%)
	B2	16 32	2.39 (59.75%) 2.74 (34.25%)	3.18 (79.5%) 4.66 (59.25%)	3.30 (82.5%) 5.86 (73.25%)	3.84 (96%) 7.07 (88.375%)
		64	3.17 (19.813%)	7.34 (45.875%)	10.87 (67.938%)	12.34 (77.125%)
		128	3.33 (10.406%)	8.87 (27.719%)	20.92 (65.375%)	24.70 (77.188%)
		4	1 (100%)	1 (100%)	1 (100%)	1 (100%)
		8	1.75 (87.5%)	1.97 (98.5%)	1.80 (90%)	1.99 (99.5%)
	B2	16 32	3.04 (76%) 4.58 (57.25%)	3.68 (92%) 6.47 (80.875%)	3.31 (82.75%) 5.97 (74.625%)	3.95 (98.75%) 7.87 (98.375%)
		64	7.00 (43.75%)	11.47 (71.688%)	10.59 (66.188%)	15.04 (94%)
		128	8.17 (25.531%)	15.86 (49.563%)	19.94 (62.312%)	30.14 (94.188%)
	Table

https://github.com/atlanmod/SparkTE_public

Acknowledgment

This paper disseminates results from the Lowcomote project, that received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skldowska-Curie grant agreement No 813884.

Experiments presented in this paper were carried out using the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER and several universities as well as other organizations (see https://www.grid5000.fr). These experiments also used the EnosLib library (see https://gitlab.inria.fr/discovery/enoslib).