Tarek Manasfi
email: tarek.manasfi@univ-amu.fr

Karine Lebaron
email: karine.lebaron@univ-amu.fr

Maxime Verlande
email: maxime.verlande@univ-amu.fr

Julien Dron
email: julien.dron@institut-ecocitoyen.fr

Carine Demelas
email: carine.demelas@univ-amu.fr

Laurent Vassalo
email: laurent.vassalo@univ-amu.fr

Gautier Revenko
email: gautier.revenko@institut-ecocitoyen.fr

Etienne Quivet

Jean-Luc Boudenne
email: luc.boudenne@univ-amu.fr

Occurrence and speciation of chlorination byproducts in marine waters and sediments of a semi-enclosed bay exposed to industrial chlorinated effluents

Keywords: Chlorination byproducts, Disinfection byproducts, Industrial chlorinated effluents, Cooling water, Marine sediments, Seawater pollution

à 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

The use of seawater in industrial cooling or heating is a common practice in many parts of the world. One of the primary operational problems of using seawater in such processes is biofouling which can result from the growth of microorganisms on surfaces where they form biofilms or the growth of macro-organisms such as clams.

Biofilms tend to stick to heat-exchange surfaces, thereby significantly reducing heattransfer coefficients. In some cases, excessive biofouling can lead to plugging of heat exchangers. There are several techniques for preventing biofouling of both types. The use of chlorine to avoid biofouling is among the most common procedures [START_REF] Khalanski | Chlorination chemistry and ecotoxicology of the marine cooling water systems[END_REF]. Chlorine is added either in the gaseous form or in the aqueous form of sodium hypochlorite solution, typically at doses of 0.5-1.5 mg/L (expressed as Cl 2 ) (Allonier et al., 1999a(Allonier et al., , 1999b;;[START_REF] Ma | Impacts of chlorination and heat shocks on growth, pigments and photosynthesis of Phaeodactylum tricornutum (Bacillariophyceae)[END_REF][START_REF] Khalanski | Chlorination chemistry and ecotoxicology of the marine cooling water systems[END_REF].

When chlorine is added to seawater, it reacts with organic and inorganic compounds present in seawater leading to the formation of a range of chlorination byproducts (CBPs) [START_REF] Boudjellaba | Chlorination by-product concentration levels in seawater and fish of an industrialised bay (Gulf of Fos, France) exposed to multiple chlorinated effluents[END_REF]. Factors such as the initial chlorine dose, temperature, pH, constitution of seawater and ultimately presence of contaminants (natural or anthropogenic) can influence these reactions leading to significant differences in the nature and levels of the different formed CBPs [START_REF] Allonier ; Allonier | Characterization of Chlorination By-products in Cooling Effluents of Coastal Nuclear Power Stations[END_REF]. The release of CBPs into the environment constitutes a concern from environmental and health standpoints. Even if the levels of CBPs formed in seawater tend to be low, their effects on the long term cannot be neglected since they are usually contained in large volumes of water over extended geographical area and for prolonged periods of time [START_REF] Boudjellaba | Chlorination by-product concentration levels in seawater and fish of an industrialised bay (Gulf of Fos, France) exposed to multiple chlorinated effluents[END_REF]. In this way, chlorine itself and its byproducts can constitute a threat to marine ecosystems [START_REF] Taylor | The effects of biological fouling control at coastal and estuarine power stations[END_REF][START_REF] Deng | Chronic exposure to environmental levels of tribromophenol impairs zebrafish reproduction[END_REF][START_REF] Pignata | Chlorination in a wastewater treatment plant: acute toxicity effects of the effluent and of the recipient water body[END_REF][START_REF] Khalanski | Chlorination chemistry and ecotoxicology of the marine cooling water systems[END_REF]. Furthermore, the spectrum of possible adverse risks is not limited to the environment since potential volatilization and subsequent photolysis of CBPs into reactive species in the atmosphere may affect human health upon exposure to CBPs or their products [START_REF] Quack | Air-sea flux of bromoform: controls, rates, and implications[END_REF]Parinet et al., 2012). Most of studies conducted previously on CBPs in the marine environment focused on nuclear and thermal power plants [START_REF] Allonier ; Allonier | Characterization of Chlorination By-products in Cooling Effluents of Coastal Nuclear Power Stations[END_REF][START_REF] Khalanski | Chlorination chemistry and ecotoxicology of the marine cooling water systems[END_REF]. These investigations looked into a limited number of compounds in cooling water effluents, often at a single or a few discharge points in the open coast. To date, data about the contamination of seawater and marine sediments exposed to multiple industrial chlorinated effluents remain very scarce. The present study investigated the contamination of seawater and sediments by CBPs in the Gulf of Fos, Southeastern France. The Gulf of Fos is a semi-enclosed bay that favors water confinement in some of its back-ends and receives the plumes of the second greatest Mediterranean river, namely Rhône river, among other freshwater inputs [START_REF] Ulses | Circulation in a semienclosed bay under influence of strong freshwater input[END_REF]. The Gulf of Fos hosts the largest port of trade in France and in the Mediterranean Sea (Marseille-Fos Port) along with a major industrial zone that includes steel, petrochemical, waste incineration, and cement industries as well as gas and electricity power plants. The aim of the study was to characterize the contamination of the bay by halogenated CBPs. For this reason, several classes of halogenated CBPs including trihalomethanes (THMs), haloacetic acids (HAAs), haloacetonitriles (HANs), trihaloacetaldehydes (THA), and halophenols (HPs) were analyzed in seawater and marine sediments obtained from the vicinity of the industrial effluents and at other sites throughout the bay. Global physicochemical parameters such as temperature, pH, salinity, bromide concentration, total organic carbon (TOC), total nitrogen (TN), and chlorine levels were also determined. Data provided by this study are crucial for any future risk assessment study for the marine pollution in the Gulf of Fos or similar marine environments.

Materials and Methods

Study Site

The study investigated the Gulf of Fos (Figure 1) located in Southeastern France, on the Mediterranean, and at about 50 km from the city of Marseille. Average water depth in the Gulf is of about 20 m. The Gulf of Fos is a semi-enclosed bay which receives several freshwater inputs with the main input being from the Rhône River.

Other minor inputs include the Berre Lagoon, irrigation and navigation canals. The In addition, there are other power plants which are irregularly active (designated by sampling stations 9x and 10x) which discharge chlorinated water (by electrochlorination). Metal industry (sampling station 11x) and oil refineries (sampling station 13x) also discharge chlorinated seawater at flows exceeding 10,000 m 3 /h.

Sample Collection

Sampling campaign was performed during spring (24, 25, and 26 April 2017). When sampling was performed, sustained southeast wind (35-40 km/h on average according to MétéoFrance, Istres) with gusts. The different sampling stations were located within the Gulf of Fos and included the main industrial outlets (sampling station names end with x) (Figure 1). Water samples were collected at all the 24 sampling stations. At each sampling station, seawater samples were collected at the surface (depth between 0 to 50 cm) and at 7m depth (or the bottom at the stations where the bottom is at depth < 7 m). Seawater samples were collected using a 5-L Niskin bottle (General Oceanics, USA). For the analysis of CBPs, sample aliquots (1 L) were placed in amber glass bottles with PTFE-lined screw caps and ascorbic acid was rapidly added to these bottles to quench any residual free chlorine. For the analysis of global physicochemical parameters including bromide ion concentration, TOC, TN, and determination of free and total chlorine, sample aliquots (1 L) were placed in amber glass bottles with PTFE-lined screw caps to which no ascorbic acid was added.

Samples were stored at 4 °C away from sunlight and extracted within 24 h from collection. Physicochemical parameters including pH, temperature, and salinity were determined on-site using a CTD-type multi-parameter probe (MS5, OTT Hydrolab, Germany).

Sediments were collected using an Ekman grab sampler at 10 sampling stations (8p, 11x, 12x, 13x, 14m, 15m, 17m, 19m, 22m, 24m). The collected sediments were stored at 4 °C away from sunlight.

Chemicals

Analytical standards including THM calibration mix, halogenated volatile mix (containing HANs), and HAA esters calibration mix, 2,3-dibromopropionic acid solution, and the THA chloral hydrate were purchased from Supelco (USA). The brominated THA tribromoacetaldehyde (97%) was purchased from Aldrich (United Kingdom) and was used to generate its hydrated form bromal hydrate (BH) in ultrapure water (Millipore, resistivity >18 MΩ cm). HPs 2-bromo-4-chlorophenol (98%) and 2,6-dibromophenol (99%) were purchased from Alfa Aesar (Germany), and 2,4-dibromophenol (95%) and 2,4,6-tribromophenol (99%) were purchased from Sigma-Aldrich (USA). A standard stock solution of each compound was prepared in methyl tert-butyl ether (MTBE, purity 99.8%) which was purchased from Merck (Germany). L-ascorbic acid, crystalline, reagent grade was purchased from Sigma (China). Sulfuric acid, analytical grade reagent, was purchased from Fisher Scientific (UK). For plotting the calibration curve, artificial seawater was spiked with the mother solutions at different concentrations and the resulting solutions were treated according to the methods described hereby for samples. Artificial seawater (ASW) was prepared according ASTM International standard practice for the preparation of substitute ocean water (method D1141-98, 2013).

Seawater sample preparation

For the analysis of THMs, HANs, THAs sample aliquots (50 mL) were first adjusted to a pH value ranging between 4.5 and 5.5 by adding sulfuric acid. For the different chemical classes, specific internal standards were added. For THMs, HANs and BH, the internal standard 1,2,3-trichloropropane was used according to U.S.EPA 551 .1 (Munch and Hautman, 1995). Samples were then extracted by liquid-liquid extraction (LLE) by adding MTBE (5 mL) and shaking manually for 3 min. Then, the organic phase was collected for analysis. For the analysis of HAAs, U.S.EPA method 552.3 [START_REF] Domino | Method 552.3 Determination of haloacetic acid and dalapon in drinking water by liquid-liquid microextraction, derivatization, and gas chromatography with electron capture detection[END_REF] was used with slight modifications. In brief, sample aliquots (40 mL) were acidified to a pH < 1 by adding concentrated sulfuric acid and extracted with MTBE (4 mL). 2,3-dibromopropionic acid was added to the extracts as a surrogate. After LLE, the organic phase containing the HAAs was collected and transferred into 15 mL vials to which acidified methanol was added and placed in a water bath at 50 °C for 2 h for derivation (methylation). The vials were then cooled, and 4 mL of saturated sodium bicarbonate solution were added before collecting the organic phase containing the HAA esters in chromatographic vials. For the analysis of HPs, derivation (acetylation) and extraction by LLE were conducted as described previously by [START_REF] Allonier ; Allonier | Characterization of Chlorination By-products in Cooling Effluents of Coastal Nuclear Power Stations[END_REF] with some modifications. In brief, 50 mL samples were mixed manually during 4 minutes with 10 g of sodium carbonate and 5 mL of acetic anhydride to derivatize the HPs. Samples were then extracted with 2.5 mL of MTBE containing the internal standard 2,4,6 trichlorophenol. The organic phases were then collected and dried with sodium sulfate before analysis.

Analytical methods

Free residual chlorine and total chlorine were measured by the colorimetric DPD method using a portable spectrophotometer (AQUALYTIC-AL 800, Germany).

Bromide levels in water were measured by an ICS-3000 Dionex ion chromatography system using a 30 mM NaOH eluent with a flow rate of 1.5 ml/min at 30 °C. Total Organic Carbon (TOC) and Total Nitrogen (TN) were measured using high temperature catalytic oxidation technique (Multi N/C 2100, Analytik Jena, Germany).

The pre-treated sample was injected (50 µl) into the furnace filled with a Pt preconditioned catalyst. The combustion was realized at 800 °C and the combustion products were carried by high purity oxygen (Linde Gas) allowing detection of CO 2 by non-dispersive infrared (NDIR) and detection of NO by chemiluminescence (CLD). Organic extracts containing CBPs were analyzed using a gas chromatograph coupled to a 63 Ni electron-capture detector (GC-ECD model Clarus 580, Perkin Elmer, Norwalk, CT, USA). An Elite 5MS capillary column was used for the separation. Helium 5.0 was used as a carrier gas at 1 mL/min. Nitrogen was used as a make-up gas at 30 mL/min. For the analysis of THMs, HANs, HPs and BH the temperature program was as follows: initially 35 °C increasing to 145 °C at a rate of 10°C/min, then at a rate of 20 °C/min up to 225 °C and finally at 10 °C/min to 260 °C, held for 2 min. For the analysis of HAAs, the temperature was initially set to 40 °C, then increased to 75 °C at a rate of 15 °C/min, then increased to 100 °C at 5 °C/min, and finally temperature reached 135 °C at 10 °C/min which and held 2 min.

Analytes were qualified using procedural standard calibration. Calibrations were performed at concentrations starting from 10 ng/L to 10 µg/L. At each concentration order, external calibrations were performed using a set of 7 standard solutions. The solutions were prepared by adding aliquots of the standard stock solution in artificial seawater and later conducting the same treatment as that of samples. Seawater and purified water reagent blanks were included with each sequence. Analysis of samples was conducted in duplicate. The detection limits (DL) and the quantification limits (QL) for the analyzed chemicals and parameters and their estimation procedure are presented in supplementary information (Table S1).

Maps were produced using the software: R (R Core Team, 2015) and Inkscape (Inkscape, 2015).

Sediment sample preparation

Sediments were brought to room temperature. Extraneous material was removed prior to homogenization. Sediment samples were homogenized in pre-cleaned collection jars by stirring vigorously with stainless steel spatulas. Dry weights were determined by placing sample aliquots in an oven at 105 °C and weighed at intervals of 24 h until successive weight differences became less than 4%. For the analysis of HPs in sediments, the method was inspired from [START_REF] Lampi | Chlorophenols in lake bottom sediments: A retrospective study of drinking water contamination[END_REF] with modifications. An aliquot (20 g) of wet sediments was weighed, 50 mL of 1M NaOH solution and internal standard (50 µL) of 2,4,6-trichlorophenol solution (1 mg/L) were added.

Then, 50 mL of hexane were added to the sample and the whole kept in an ultrasonic bath for 10 min. The hexane phase was then discarded, and the aqueous phase was collected and introduced into 65-mL glass vials with PTFE-lined screw caps which were centrifuged at 2000 rpm for 10 min. 40 mL of the resulting supernatant was then collected and was treated according to protocols used for the analysis of HPs in water samples described above.

Results and Discussion

Physicochemical parameters

At the different sampling stations, physicochemical parameters including water temperature, pH, salinity, TOC, TN, bromide content, and levels of free residual chlorine, and total chlorine were determined. Table S2 presents the measurements of the physicochemical parameters at the sampling stations. Very little variation was observed in the pH among the different sampling points located throughout the gulf (pH ranged between 8.13 to 8.25). stations especially at the surface freshwater intrusion appeared at stations lower salinities than the sampling stations bottom (or at 7 m depth) the salinity was generally homogeneous (> 37 remained slightly lower in the eastern sector. located on the canal receiving water of the Rhône river, had unusually high salinities (24.2 and 34.2 g/L on the surface and bottom, respectively) resulting from the vigorous mixing of the Rhône coming from the Gulf of Fos especially that on the day of sampling moderate to strong southeast winds were registered (pH ranged between 8.13 to 8.25). Salinity varied among the different sampling surface (Figure 2). A marked influence of the Rhone River appeared at stations located nearby (2c, 8p and 5p lower salinities than the sampling stations located offshore.

Variation of surface water salinity (in g/L) at the different sampling stations seems to drive these less salty waters to the western part of the Gulf as appeared at the surface at stations 17m, 18m, 19m, 21m and 22m m depth) the salinity was generally homogeneous (> 37 slightly lower in the eastern sector. The sampling station 2 canal receiving water of the Rhône river, had unusually high salinities (24.2 and 34.2 g/L on the surface and bottom, respectively). This salinity is probably resulting from the vigorous mixing of the Rhône River freshwater with seawater ulf of Fos especially that on the day of sampling moderate to strong southeast winds were registered. Bromide concentrations were generally among the different sampling marked influence of the Rhone River 8p and 5p) which had at the different sampling stations less salty waters to the western part of the Gulf as appeared at the surface at stations 17m, 18m, 19m, 21m and 22m. At the m depth) the salinity was generally homogeneous (> 37 g/L), but he sampling station 2c, which was canal receiving water of the Rhône river, had unusually high salinities . This salinity is probably freshwater with seawater ulf of Fos especially that on the day of sampling moderate to Bromide concentrations were generally homogenous among the different sampling stations within the Gulf (Table S2).

Bromide levels were relatively high even at points close to the intrusion of freshwater of the Rhone River. Bromide levels at sampling stations 3p to 8p ranged from 46.6 to 58.21 at the surface, and from 53.75 to 69.34 at the bottom/depth. Temperature at the surface and in deep water are presented in Table S2. Excluding outlets, water temperatures at the bottom were globally homogeneous across the Gulf. Water temperature was 14.5 to 15 ° C at the surface and 14 ° C at the bottom/depth. This very little gradient is compatible with the agitated conditions that favor the mixing of water bodies. The temperature of freshwater (1c, 2c) was slightly higher (15 to 16 °C). Near the outlets at sampling stations 9x, 10x, 11x, temperatures were slightly higher at the surface (up to 18.9 °C) and the bottom (up to 19.8 °C), due to the discharge of heated seawater from the nearby industries. In general, the surface temperatures at the sampling stations close to the shore were higher than offshore stations (Figure 3). The levels of the CBPs identified in the Gulf of Fos are presented in Table S3. Several CBPs belonging to different chemical classes were detected in the Gulf. The main detected CBPs included bromoform, dibromoacetic acid (DBAA), dibromoacetonitrile (DBAN), tribromoacetic acid (TBAA), and BH. Low concentrations of bromochloroacetonitrile (BCAN) were also detected near some outlets (sampling stations 8p, 9x, and 10x). Chloroform was detected at some sampling stations at levels below the QL. The occurrence of predominantly brominated CBPs in chlorinated seawater can be explained by the formation of bromine upon the addition of chlorine to seawater. This speciation is in agreement with previous studies that reported the formation of brominated byproducts in chlorinated bromide-containing water [START_REF] Ged | Effect of seawater intrusion on formation of brominecontaining trihalomethanes and haloacetic acids during chlorination[END_REF]. In the presence of appreciable amounts of bromide ions, chlorine oxidizes bromide ions and forms hypobromous acid and hypobromite ions [START_REF] Singer | Humic substances as precursors for potentially harmful disinfection by-products[END_REF]. As an oxidant, bromine is stronger than chlorine and reacts 10 times faster with organic matter (Westerhoff et al., 2004). The reactions leading to the formation of bromine in bromide-rich water such as seawater are as follows [START_REF] Heeb | Oxidative treatment of bromide-containing waters: Formation of bromine and its reactions with inorganic and organic compounds -A critical review[END_REF]:

HOCl + Br -→ HOBr + Cl -k 1 = (1.55-6.84).10 3 M -1 s -1 (1)

ClO -+ Br -→ BrO -+ Cl - k 2 = 9.10 -4 M -1 s -1 (2)
Although HOBr and OBr -are the most abundant species involved in these reactions, several bromine species can react with organic compounds present in seawater [START_REF] Heeb | Oxidative treatment of bromide-containing waters: Formation of bromine and its reactions with inorganic and organic compounds -A critical review[END_REF]Manasfi et al., 2017). Bromine species such as Br 2 , Br 2 O, BrOCl, and BrCl are less abundant but more reactive than HOBr and OBr -and have been shown

to contribute to bromination of some organic compounds [START_REF] Sivey | Reactivity of BrCl, Br2, BrOCl, Br2O, and HOBr toward dimethenamid in solutions of bromide + aqueous free chlorine[END_REF]. Table 1 shows the equilibrium concentrations of bromine species when seawater containing 67 mg/L bromide, at pH 8.05, is chlorinated with active chlorine (2 mg/L). These concentrations were estimated using PHREEQC [START_REF] Parkhurst | Description of input and examples for PHREEQC version 3: a computer program for speciation, batch-reaction, onedimensional transport, and inverse geochemical calculations[END_REF].

The levels of CBPs were highest at the sampling stations located near the outlets where chlorinated industrial effluents are released (stations 8p, 9x, 10x, 11x, and 13x). At these outlets, the concentration of bromoform ranged from 0.17 to 1.95 µg/L (at the surface) and from 0.46 to 2.36 µg/L (at the bottom). At the sampling station 12x (in the vicinity of the outlet of LNG Fos-Cavaou terminal), surprisingly low concentrations were detected. This aberration is probably due to the dislocation of the sampling Niskin bottle under the effect of wind during sampling. At its exit at sea, station 12x, at 150 m from the shore, strong wind and current of south-east, made the measurements and sampling very difficult. Throughout the studied area of the Gulf, bromoform concentrations ranged from 0.05 to 1.95 µg/L and from 0.06 to 2.3 µg/L at the surface and the bottom, respectively. The slightly higher concentrations of bromoform at the bottom compared to the surface may result from accelerated volatilization on the surface because of wind and slightly higher temperatures. For the other CBPs, distinctive discrepancies between bottom and surface levels were not observed. DBAA concentrations ranged from 0.35 to 1.40 µg/L and from 0.39 to 1.23 µg/L at the surface and the bottom, respectively. DBAN concentrations ranged from 0.049 to 0.67 µg/L and from 0.05 to 0.79 µg/L at the surface and the bottom, respectively. BCAN was detected at concentrations ranging from 0.04 to 0.05 µg/L at the surface and the bottom. BH concentrations ranged from 0.075 to 0.13 µg/L and from 0.06 to 0.18 µg/L at the surface and the bottom, respectively. Overall, the concentrations of CBPs detected in the Gulf of Fos in the present study were relatively low compared to levels previously reported for these compounds in seawater exposed to industrial chlorinated effluents. In a survey conducted in coastal power stations, at the point of discharge, bromoform concentrations varied from 3.1 µg/L to 29.20 µg/L and DBAN concentrations varied from 0.10 µg/L to 11.39 [START_REF] Khalanski | Chlorination chemistry and ecotoxicology of the marine cooling water systems[END_REF]. Since many of the CBPs are volatile, increased volatilization with accentuated winds could be one reason. Another reason could be related to decrease in chlorination activity at some of the plants since water temperature were not very high which naturally decreases biofouling.

Despite the low concentrations, bromoform was still detected in the Gulf of Fos even at sampling stations that were relatively distant from the outlets (such as 14m, 15m, 16m, 17m, 19m, 20m, and 24m). The concentrations of bromoform at these sampling stations away from the outlets ranged from 0.06 to 0.57 µg/L and from 0.06 to 0.64 µg/L at the surface and the bottom respectively. These concentrations are superior to the typical background levels of bromoform emitted by marine algae in seawater that is not exposed to chlorinated effluents. The latter has been estimated at 0.025 μg/L, and rarely exceed 0.1 μg/L unless when extensive beds of macro-algae are present, which is not the case in the Gulf of Fos [START_REF] Quack | Air-sea flux of bromoform: controls, rates, and implications[END_REF]. The diffusion of bromoform across the Gulf of Fos is represented in Figure 4. This figure shows that shows contamination of the Gulf of Fos with bromoform at concentrations lower than what is detected near the outlets, yet higher than what is found in uncontaminated seawater.

Distribution of CBPs in seawater

Among the detected CBPs, bromoform was the most abundant species at most of the sampling stations (Figure 5).

This finding is in agreement with previous studies that investigated the occurrence of CBPs in seawaters exposed to chlorinated industrial effluents [START_REF] Allonier ; Allonier | Characterization of Chlorination By-products in Cooling Effluents of Coastal Nuclear Power Stations[END_REF][START_REF] Boudjellaba | Chlorination by-product concentration levels in seawater and fish of an industrialised bay (Gulf of Fos, France) exposed to multiple chlorinated effluents[END_REF][START_REF] Jenner | Chlorination byproducts in chlorinated cooling water of some European coastal power stations[END_REF].

Another THM, dibromochloromethane was detected at levels below its QL. DBAA was the second most abundant CBP at all the sampling sites except station 11x where DBAN constituted the second most abundant species. This discrepancy can be attributed to the presence of high levels of nitrogen-containing organic compounds near the industrial outlet 11x. Furthermore, pH is another factor that has been shown to affect the speciation of CBPs [START_REF] Cimetiere | Effect of some parameters on the formation of chloroform during chloramination of aqueous solutions of resorcinol[END_REF]. [START_REF] Hansen | Effect of pH on the formation of disinfection byproducts in swimming pool water -Is less THM better?[END_REF] demonstrated that with decreasing pH, THM formation was reduced while HAN formation increased and HAA formation remained constant. However, in this study the influence of pH on the distribution of CBPs can be ruled out since the pH was largely homogenous across the Gulf (pH = 8.1 -8.2). While DBAA was generally the second most abundant CBP, the tribrominated HAA TBAA was detected only at the sampling station 8p. This finding may be explained by the lower stability of TBAA in comparison to DBAA [START_REF] Manasfi | Identification of disinfection by-products in freshwater and seawater swimming pools and evaluation of genotoxicity[END_REF][START_REF] Zhang | Decomposition of trihaloacetic acids and formation of the corresponding trihalomethanes in drinking water[END_REF]. It has been reported that TBAA may decompose to form bromoform in aqueous solutions [START_REF] Zhang | Decomposition of trihaloacetic acids and formation of the corresponding trihalomethanes in drinking water[END_REF]. HPs were not detected in seawater samples.

Occurrence of brominated halophenols in sediments

Sediments were collected at the sampling stations (8p, 11x, 12x, 13x, 14m, 15m, 17m, 19m, 22m, 24m). Two bromine-containing HPs were detected at the sampling stations 22m and 13x at concentrations in the order of ng/g level (of dry weight). At sampling station 22m, 2-chloro-4-bromophenol and 2,4,6-tribromophenol were detected at 1.8 and 2.1 ng/g (of dry weight). At sampling station 13x, the two HPs were detected at 0.3 and 1.5 ng/g (dry weight), respectively. To our knowledge, this is the first study to report the levels of some halophenols in marine sediments obtained from the vicinity of industrial chlorinated effluents. Although the detected HPs were brominated and therefore compatible with the general speciation of CBPs determined in seawater in the Gulf, further investigations are necessary to discriminate whether the detected halophenols originate from the chlorination of seawater or are produced from algal species and biota.

Conclusions

Seawater in the Gulf of Fos which is exposed to multiple industrial chlorinated

Conflicts of interest

Authors declare no conflicts of interests.

Fig. 1 .

 1 Fig. 1. Overview of the Gulf of Fos 24m).

Figure 2 .

 2 Figure 2. Variation of surface water salinity (

Figure 3 .

 3 Figure 3. Variation of surface water temperature (in °C) at the different sampling stations

  effluents contained predominantly the brominated CBPs bromoform, DBAA, DBAN, and BH. Other CBPs including TBAA and BCAN were also detected at some sampling stations. The levels of CBPs were relatively low compared to levels reported in the literature for seawater exposed to chlorinated effluents from nuclear plants and other industries which employ chlorination. Despite the low levels, bromoform was found at sampling stations relatively far from the chlorinated effluent outlets at concentrations superior to what is found naturally in uncontaminated seawater. In sediments, two HPs namely 2-chloro-4-bromophenol and 2,4,6-tribromophenol were detected at two sampling stations. Since the concentrations of CBPs depend on many factors including meteorological conditions, season, and chlorination operations in the nearby industrial plants, future investigations that assess the levels of CBPs under different meteorological conditions and season are warranted. Furthermore, the importance of risk assessment studies that estimate the impact of the detected concentrations of CBPs on the environment is highlighted.

Acknowledgements

This work was funded by the French National Research Agency(ANR) within the project FOSSEA (ANR-16-CE34-0009). Authors acknowledge Météo-France for meteorological data through the agreement N° DIRSE/REC/16/02/0 and reviewers for helpful comments.


Kaia Tõnsuaadu
email: kaia.tonsuaadu@ttu.ee

Michel Gruselle

Frieda Kriisa

Andres Trikkel

Patrick Gredin

Didier Villemin

Dependence of the interaction mechanisms between L-serine and O-phospho-L-serine with calcium hydroxyapatite and copper modified hydroxyapatite in relation with the acidity of aqueous medium

Keywords: HAp, hydroxyapatite, CaHAp, calcium-hydroxyapatite, CuHAp, copper-substituted hydroxyapatite, L-Ser, L-serine, Calcium hydroxyapatite, Copper modified hydroxyapatite, L-Serine, O-phospho-L-serine, adsorption

come

Dependence of the interaction mechanisms between l-serine and O-phospho-l-serine with calcium hydroxyapatite and copper modified hydroxyapatite in relation with the acidity of aqueous medium

Introduction

Hybrid organic-inorganic materials are of first importance in biological processes. Due to the role of apatites in the growth of human bones and teeth and in drug delivery, there is an increasing interest in studies related to the interactions between amino-acids and calcium apatites. Many works are devoted to the adsorption of organic and bioorganic molecules on the surface of apatites [START_REF] Šupová | Substituted hydroxyapatites for biomedical applications: A review[END_REF][START_REF] El Rhilassi | Interaction of some essential amino acids with synthesized poorly crystalline hydroxyapatite[END_REF][START_REF] Errassifi | Advances in bioceramics and biotechnologies[END_REF][START_REF] Kollath | A modular approach to study protein adsorption on surface modified hydroxyapatite[END_REF][START_REF] Errassifi | Infrared, Raman and NMR investigations of risedronate adsorption on nanocrystalline apatites[END_REF]. In general, these experiments, carried out at pH higher than 6, point out the fact that apatite surface is not fundamentally modified during the adsorption process, specifically, release of mineral ions like calcium or phosphate ions from the apatite structure are not taken into consideration in the interactions. Therefore, the role of bioorganic molecules, in particular, amino-acids (AA) as reagents to change composition of the mineral support, is neglected.

Particular interest among other AA has been paid on L-Serine (L-Ser) and O-phospho-L-serine (O-Ph-L-Ser) adsorption, which were characterized by Langmuir-type adsorption on apatite [START_REF] Spanos | Model Studies on the Interaction of Amino Acids with Biominerals: The Effect of L-Serine at the Hydroxyapatite-Water Interface[END_REF][START_REF] Benaziz | Adsorption of O-Phospho-L-Serine and L-Serine onto Poorly Crystalline Apatite[END_REF][START_REF] Skartsila | Adsorption monitoring of phospho-l-serine on hydroxyapatite[END_REF][START_REF] Spanos | Model Studies of the Effect of Orthophospho-L-serine on Biological Mineralization[END_REF][START_REF] Misra | Interaction of ortho-Phospho-l-serine with Hydroxyapatite: Formation of a Surface Complex[END_REF]. The studies carried out at pH 7-10 with L-Ser led to the conclusion that L-Ser adsorbs on the surface of hydroxyapatite (HAp) through electrostatic attractions exerted between one negative site of the HAp surface, i.e., phosphate or hydroxyl ion, and the positively charged protonated amino group of one Ser molecule, forming an ion pair surface complex [START_REF] Spanos | Model Studies on the Interaction of Amino Acids with Biominerals: The Effect of L-Serine at the Hydroxyapatite-Water Interface[END_REF][START_REF] Benaziz | Adsorption of O-Phospho-L-Serine and L-Serine onto Poorly Crystalline Apatite[END_REF].

O-Ph-L-Ser has higher affinity for the HAp surface compared to L-Ser molecules [START_REF] Benaziz | Adsorption of O-Phospho-L-Serine and L-Serine onto Poorly Crystalline Apatite[END_REF]. The enhanced adsorption capacity noted for O-Ph-L-Ser might be related to the presence of phosphate groups in the molecule, which are potential specific attachment sites. The monoprotonated negatively charged (HL 2-) O-Ph-L-Ser species were found adsorbed forming surface complexes with the positively charged ≡CaOH 2+ sites on the surface of HAp at pH 7 [START_REF] Skartsila | Adsorption monitoring of phospho-l-serine on hydroxyapatite[END_REF]. The negatively charged deprotonated carboxyl and phosphate groups of the adsorbed HL 2-species orient themselves at the maximum possible distance from the negatively charged surface of HAp because of electrostatic repulsions [START_REF] Spanos | Model Studies of the Effect of Orthophospho-L-serine on Biological Mineralization[END_REF].

The only study performed in acidic O-Ph-L-Ser solutions [START_REF] Misra | Interaction of ortho-Phospho-l-serine with Hydroxyapatite: Formation of a Surface Complex[END_REF] concluded that there occurs an ion-exchange between HAp and P-Ser in dilute solutions. The dissolution of apatite released phosphate and calcium ions into the solution in the amounts linearly related to the uptake of serine (<8 mM solution).

Taking into consideration that the drug delivery process is strongly influenced by the acidity of the biological medium, we decided to revisit the mechanisms of the adsorption of AA, namely L-Serine and O-Phospho-L-serine, shown in Fig. 1, with Ca-hydroxyapatite (CaHAp) and a copper modified Ca-hydroxyapatite (CuHAp) in relation with the acidity of the aqueous medium from acidic to smoothly basic one. The use of calcium apatites doped with copper ions is interesting to show the specific role that oxidation state II ions can play in sorption mechanism of AA. The question is therefore whether these ions will increase the ability of apatite to adsorb amino acids on the apatitic surface and by what mechanism. We first chose copper ions for their known role in biological and catalytic processes [START_REF] Festa | Copper: An essential metal in biology[END_REF][START_REF] Clavadetscher | Copper Catalysis in Living Systems and In Situ Drug Synthesis[END_REF][START_REF] Chemler | Copper catalysis in organic synthesis[END_REF].

Experimental Section

Synthesis of CaHAp and CuHAp

The chemicals used for synthesis of apatites were of ACS grade, purchased from Sigma-Aldrich, Merck KGaA and Acros. The CaHAp was synthesized by precipitation in aqueous solution using a PC-controlled LARA reactor system (Radleys) in deionized CO 2 -free water at 25 °C in N 2 atmosphere keeping pH=9 by continuous addition of ammonia solution and keeping Ca/P mole ratio equal to 1.67. The Ca(NO 3 ) 2 (0.5 M) and NH 4 H 2 PO 4 (0.3 M) solutions were added simultaneously into NH 4 OH aqueous solution (0.5 L). The obtained suspension was stirred for 1 hour and matured for 24 hours at room temperature. The precipitate was filtered, washed several times with deionized CO 2 -free water and lyophilized. The Cu substituted apatite was obtained mixing the washed CaHAp slurry in Cu-acetate solution (0.01 M) for 3 hours, after that filtrated, washed and lyophilized. Apatite, 300 mg, was added to 50 ml of AA solution, mixed in centrifugal tube at overhead mixer for 100 min and centrifuged. The experiments were performed in triplicate and the results were averaged. After separation of the solid and liquid phases, pH and the cations content were determined in the solution and the UV spectrum was recorded. The solid residue was dried at 60 °C and then subjected to XRD, FTIR and thermal analysis.

Analytical methods

XRD patterns were recorded with an INEL Equinox 1000 diffractometer operating in real time acquisition over 2θ range 5 -115° in asymmetric mode. The diffractometer is equipped with a cobalt anode and a primary graphite monochromator ( λ Kα1 = 1.7889 Å, λ Kα2 = 1.7928 Å). The 2θ angular resolution is 0.029. The cell dimensions and average crystallite size were determined by FULLPROF program [START_REF] Rodriguez-Carvajal | Resent developments of the program FULLPROF[END_REF].

FTIR spectra were recorded with a Bruker 55/S/NIR FTIR ATR spectrometer as an average of 30 scans at resolution 4 cm -1 .

The simultaneous TG/DTA/EGA-MS measurements were performed in an apparatus consisting of a LabSys-Evo 1600 (Setaram, France) thermal analyzer and an OmniStar (Pfeiffer, Germany) Quadrupole mass spectrometer. Coupling between the two components was provided through a heated capillary tube kept at T = 180 °C. The ion currents of the selected mass/charge (m/z) numbers were monitored in multiple ion detection (MID) mode (Quadera version 4.20 software) with the collection time of 1 s for each channel. The measurement was carried out in Ar-O 2 (20%) mixture atmosphere in the temperature range of 30-1000°C using the heating rate of 10 deg•min -1 , the gas flow rate of 60 mL•min -1 and open Pt crucibles. The initial mass of sample was 25-35 mg. The MS mass-to-charge ratios (m/z) selected for analysis were 15 for NH 3 to avoid the overlapping with H 2 O, 18 for H 2 O, 30 for NO x (sum of NO, N 2 O and NO 2 ) and 44 amu for CO 2 .

UV absorption spectra of the solutions in deionized water were recorded with Biochrom Libras 70PC UV/Visible Spectrophotometer in the wavelength interval from 190 to 400 nm with step 2 nm in quartz cuvettes of 10 mm path length using deionized water as reference.

The macro-components of the synthesis products were determined by standard analytical methods. Content of Ca and Cu was determined after dissolution of sample in diluted HCl by The error in chemical analyses was under 2 %.

The specific surface area (SSA) measurements were performed by BET-method (adsorptive gas N 2 , carrier gas He, heating temperature 110°C) using sorptometer EMS-53.

Results

Characterization of apatites

The XRD, FTIR and chemical analysis revealed that the materials obtained by precipitation were near-stoichiometric apatites, where cations/P mole ratio was 1.62 (Table 1, theory 1.67).

XRD patterns and the unit cell parameters of the precipitated apatites, CaHAp and CuHAp, correspond well to a standard HAp (Diffraction standards ICDD card 00-055-0592) [START_REF] Rodriguez-Carvajal | Resent developments of the program FULLPROF[END_REF]. The unit cell parameters were a = 9.413 and 9.394, c = 6.860 and 6.855 Å (±0.001Å), respectively.

Decrease in unit cell parameters of CuHAp is in correlation with the smaller atomic radii of Cu in comparison with Ca. The average size of particles 110 (±21) and 68 (±17) nm correlates with the specific surface values 126 and 154 m 2 /g. To analyze the consequences of the interactions between the solid material and the AA solutions, the changes in the solutions pH, the quantities of Ca and Cu ions released in the solution during contact and the chemical and structural transformations of the apatite materials were determined. The pH value changes in correlation with the AA concentration and initial pH value. The largest increase takes place at higher concentration and at initial pH <7. The pH change is almost at the same level for both AA and apatites. At initial pH 8, the decrease up to -0.9 is followed, which is more evident in L-Ser solutions.

UV spectroscopic analysis

Frequently, UV absorption measurements have been used for amino acid content determination in a solution [START_REF] Misra | Interaction of ortho-Phospho-l-serine with Hydroxyapatite: Formation of a Surface Complex[END_REF]. Unfortunately, the shape of L-Ser and O-Ph-L-Ser absorption spectra depends not only on AA concentration but also on the content of other ions (NO 3 ̄, PO 4 3¯̄, Cu 2+ ) and pH of the solution (see Suppl. Mat. Fig. S3-S6). Therefore, the determination of amino acid concentration by UV spectra is not reliable in the conditions when several variables as pH and ions concentration are changing.

However, UV spectra revealed the formation of Cu-AA complexes in the solutions. In Fig. 2.

UV spectra of L-Ser containing 0.0015 mg Cu per 1 ml at different pH values (A) and the intensities of the absorption peak at 246 nm depending on pH value of Cu-acetate, L-Ser + Cuac., and O-Ph-L-Ser + Cu-ac. solutions (B) are presented. The intensity of the peak at 246 nm that corresponds to Cu-AA complex, increases equally for both AA up to pH 6 and keeps this value at higher pH. The peak intensity at 246 nm depends also on Cu 2+ content in solution (Suppl. Mat. Fig. S3). Therefore, Cu 2+ release from apatite is related to its complexation with AA. S1.

The change of unit cell parameters as a result of interaction with L-Ser and O-Ph-L-Ser solutions indicates the small changes in apatite crystal structure. Due to the reaction with L-Ser, both unit cell parameters a and c increase while cell parameters of apatite remain almost unchanged or decrease very slightly by interaction with O-Ph-L-Ser. Actually, these changes are too small to be correctly interpreted.

FTIR spectroscopy

Comparison of the spectra recorded before and after amino acid sorption experiments reveals very small changes in the apatite spectra. These spectra are shown in Figures 3 and4. However, some differences are found at wavelength intervals of 700-850 and 1300-1700 cm -1 .

In the spectrum of a mechanical mixture of CaHAp and L-Ser (5 %) the characteristic absorptions of L-Ser and apatite are found at the same wavelength values that in the pure forms: at 1586 cm -1 for CO 2 -(ν as ), at 1417 cm -1 for CO 2 -(ν s ), at 1600 cm -1 for COO -(ν s ) , at 1464 cm -1

for NH 3 + (ν s ), and at 1306 and 1342 cm -1 for CH groups vibrations [START_REF] Jarmelo | Infrared and Raman spectroscopic characterization of the hydrogen-bonding network in L-serine crystal[END_REF]. In the spectra of the material resulting from the sorption of L-Ser on apatite surface, an absorption peak at 1547 cm -1 that could be assigned to COO -(ν as ) [START_REF] Lee | High protein adsorptive capacity of amino acid-functionalized hydroxyapatite[END_REF] is visible. The absorption of carbonate substituted for phosphate in apatite structure [START_REF] Fleet | Carbonated Hydroxyapatite: Materials, Synthesis, and Applications[END_REF] or of NH 3 + (ν s ) at 1453 cm -1 [START_REF] Benaziz | Adsorption of O-Phospho-L-Serine and L-Serine onto Poorly Crystalline Apatite[END_REF] is increased. A shoulder at 1630 cm -1 appears next to H 2 O absorption at 1641 cm -1 in apatite structure, which has been assigned to COO -(ν as ) vibration shifted due to L-Ser adsorption on apatite [START_REF] Gonzalez-Mcquire | Synthesis and characterization of amino acid-functionalized hydroxyapatite nanorods[END_REF] and absorption at 1329 cm -1 appears in spectra of samples treated at pH 2. These changes can be seen better for the samples of CaHAp. No impact of the solution concentration is found. particularly after reaction at pH 2, and the increase in peak intensity at 1420 cm -1 , related to COO -(ν s ) vibration [START_REF] Benaziz | Adsorption of O-Phospho-L-Serine and L-Serine onto Poorly Crystalline Apatite[END_REF][START_REF] Lee | High protein adsorptive capacity of amino acid-functionalized hydroxyapatite[END_REF], are followed in the spectra of O-Ph-L-Ser treated apatites.

Additionally, to the changes in the spectral region of 1300-1700 cm -1 , the peak at 745 cm -1 , assigned to P-C (ν s ) [START_REF] Othmani | Synthesis and characterization of hydroxyapatite-based nanocomposites by the functionalization of hydroxyapatite nanoparticles with phosphonic acids[END_REF], is shifted up to 762 cm -1 depending on the final pH of the solution.

This peak intensity increases together with the O-Ph-L-Ser concentration and decrease in the initial pH increase of the solution, similarly, to Ca release. The absorption peaks related to NH 3 + vibrations are less intensive but the changes related to COO -vibration intensities are more visible.

Thermal analysis

The thermal analysis (thermogravimetric analysis together with thermal changes and evolved gas analyses) were performed with the aim to prove the adsorption of AA on apatite and also to evaluate the amount of AA bound. The amount of organic matter in a solid sample was determined from the amount of NO x gases evolved at calcination detected by MS analysis. The peak area under the NO x evolvement curve was measured and the normalized value was used as an indicator of the amino acid content in the sample.

The mechanical mixtures of apatite and O-Ph-L-Ser or L-Ser with 3 and 5 mass % in it were prepared for estimating the NO x amount released at calcination from sorption experiment samples. The results of TA are given in Figures 567. TA clearly exhibits the adsorption of L-Ser and O-Ph-L-Ser on apatite by an exothermic effect with the maximum in the interval of 200-300 °C caused by the oxidation reaction of organic matter. In addition, the mass loss and the amount of CO 2 and NO x released, increases in comparison with the pure apatites (Fig. 5). The profile of the curves changes together with the variation in the solution pH used (Fig. 6), indicating to different sorption mechanism at different pH. The amount of NO x evolved increases with the increase in AA solution concentration and decreases with the increase in pH (Fig. 7). By TA results the amount of AA adsorbed on CuHAp is almost the same as on CaHAp.

Discussion

To clarify phenomena at the apatite surface when opposed to AA at different pH in aqueous medium, the behaviour of both components must be taken into consideration:

1) The apatite surface properties, which depend on cations/P mole ratio and chemical composition and on the medium pH [START_REF] Tsuchida | Reaction of ethanol over hydroxyapatite affected by Ca/P ratio of catalyst[END_REF][START_REF] Bengtsson | Surface complexation and proton-promoted dissolution in aqueous apatite systems[END_REF][START_REF] Rodenas | Surface properties of various powdered hydroxyapatites[END_REF].

2) The ionization state of AA, which is strongly dependent on the medium acidity [START_REF] Fleck | Salts of amino acids: crystallization, structure and properties[END_REF].

It is found that for the initial pH range 4-10, the final pH values of apatite water suspensions are close to the pH PZC , which is the pH value where the surface charge is equal to zero, namely, the pH at which the charge due to the positive surface groups is equal to negative ones [27]. For the CaHAp and CuHAp used in the study pH PZC were found to be 6.46 and 6.82, respectively.

Above the pH PZC the apatite surface is negatively charged, it is lightly positive at pH 6 and highly positive at pH 2, at which the apatite dissolves easily [START_REF] Dorozhkin | Dissolution mechanism of calcium apatites in acids: A review of literature[END_REF].

L-Ser appears in water solution at pH 2 partly in the form of neutral zwitterion and partly in the form of positively charged cation; at pH 6 as a neutral zwitterion and at pH 8 about 10 % of protons are neutralized and negatively charged -COO -groups are formed according to Scheme 1.

Scheme 1. Ser dissociation scheme [START_REF] Fleck | Salts of amino acids: crystallization, structure and properties[END_REF].

Therefore, as a result of apatite dissolution in acidic medium the pH of L-Ser solution increases and the rest of HAp [+] is in equilibrium with L-Ser zwitterion, which enables complexation via COO -groups. This explains the COO -vibrations shifts in FTIR spectra (Fig. 3).

Almost no electrostatic tension takes place between HAp and L-Ser in a solution at pH 6-8. L-Ser can adsorb on hydroxyapatite surface only due to the complexation with surface cations. The complexation constant logK f of L-Ser is 1.43 and 7.90 for Ca 2+ and Cu 2+ , respectively [START_REF]Standard Reference Database 46, NIST Critically Selected Stability Constants of Metal Complexes Database Version[END_REF]. The higher logK f of Cu 2+ explains the much higher Cu release in comparison with Ca in sorption experiments.

The changes in HAp IR spectra (appearance of a peak at 1546 cm -1 , Fig. 3), related to the COO -group complexation with Ca 2+ sites at HAp surface, in contrast with the unchanged IR spectra of CuHAp, as well as the high level of Cu 2+ ions release, support the hypothesis of the formation of soluble L-Ser-Cu complex.

The higher organic matter content in the samples after sorption experiments at lower pH values (Fig. 7) can be explained by higher positive surface charge of apatite and, accordingly, higher electrostatic tension in this case. negatively charged zwitterion electric tension exists, that results in the highest content of organic matter in the samples after sorption experiment at initial pH 2, detected by TA (Fig. 7) and more visible changes in IR spectra (Fig. 3 and4).

The electric tension and the stability of Ca-P-Ser soluble complexes decrease at higher pH values, when the stability of Cu complexes increases [START_REF] Śmiechowski | Theoretical pKa prediction of O-phosphoserine in aqueous solution[END_REF][START_REF] Zachariou | Potentiometric Investigations into the Acid-Base and Metal Ion Binding Properties of Immobilized Metal Ion Affinity Chromatographic (IMAC) Adsorbents[END_REF]. The increased stability of O-Ph-L-Ser -Cu complexes, in turn, explains the increase in Cu release at pH 8.

Therefore, the cations release from apatite takes place as a consequence of apatite dissolution at pH <5 and also due to the formation of soluble complexes with amino acids. At higher pH values the main reason of the cations release is the formation of complexes, particularly with Cu 2+ ions. In acidic medium and in the case of O-Ph-L-Ser, the effective site for bonding of metal ions is the phosphate group, which significance decreases with increasing pH as the effectiveness of the carboxyl and amine groups is raised [START_REF] Jastrzab | Effectiveness of Phosphate Groups in Noncovalent Interactions in Binary Adenosine Nucleotides/Phosphoserine Aqueous Systems[END_REF].

Conclusions

The processes occurring on the apatite surface when treated by L-Ser and O-Ph-L-Ser aqueous solutions were clarified. To understand the adsorption phenomena, it is necessary to take into consideration the ionic state of the amino acid and that of the apatite surface which are tightly correlated with the pH value of the solutions. In acidic solutions the dominant process is apatite dissolution. At the same time, due to the electrostatic attraction, negatively charged ions of the amino acid are adsorbed on the apatite surface. The amino acids interaction with apatite at pH > 6 is mainly related to its complexation with cations of apatite structure. In accordance with the higher stability of soluble Cu 2+ complexes, the relative release of Cu 2+ ions from Cu substituted apatite is higher than that of Ca 2+ ions in both amino acid solutions. Thereby, if the Cu release increases with the pH and O-Ph-L-Ser concentration increase, the Ca release drops remarkably at pH 8 in the case of P-Ser amount 1 mmol per 300 mg HAp. Impact of pH is less important for Cu release in L-Ser solution, the effect is more visible for Ca release in both AA solutions.

Fig. 1 .

 1 Fig. 1. L-Serine (left); O-Phospho-L-serine (right).

2. 2 .

 2 Sorption experiments L-Ser and O-Ph-L-Ser sorption experiments were performed with AA solutions in the range of 0.005 to 0.05 M in the pH interval from 2 to 8, adjusted by addition of 0.01 M KOH or 0.01 M HNO 3 solution using automate titrator Mettler Toledo T90 with pH electrode DG-112 Pro. The accuracy was ±0.02. The initial pH values were chosen according to the O-Ph-L-Ser solution titration curve.

  SpectrAA 55B Flame AAS) and phosphorus by photo colorimetry method as phosphomolybdate yellow complex (Biochrom Libras 70PC UV/Visible Spectrophotometer).

3. 2 .

 2 Analysis of solutions 3.2.1. Change of the AA solutions pH after contact with CaHAp and CuHAp Due to various processes that occur at the solid surface/solution interface (preferential dissolution of certain constituents of crystal lattice, ionization of surface groups, adsorption of ions or formation of complex compounds between surface groups and ions from the solution), the final pH values differ from the initial ones. After mixing HAp -water suspension during 100 min the final pH values are found as 6.46 and 6.82 for CaHAp and CuHAp, respectively. The changes of pH values of the L-Ser and O-Ph-L-Ser solutions as a result of the contact with apatite are presented in Suppl. mat. Fig. S1. The bigger change occurs at low initial pH. The pH variation goes from the initial pH value 2 to a value greater than 4. For L-Ser starting at pH 6 after mixing, the pH values increase in all the studied cases by 0.5 to 1.25 pH units. When the initial pH is 8, after mixing, the change in all pH values is -0.5 to -1.0. These trends show that the solution tends to return to the pH corresponding to the point of zero charge (pzc). For O-Ph-L-Ser the same trends are observed. The bigger change occurs at low pH, the pH variation goes from the initial pH=2 to a value greater than 4.

Fig. 2 . 3 . 3 .

 233 Fig. 2. A: UV spectra of 0.01 M L-Ser + 0.0015 mg Cu/ml at different pH. B: Absorption

Fig. 3 .

 3 Fig. 3. FTIR spectra of CaHAp and CuHAp before and after sorption experiment in 0.02 M L-

Fig. 4 .

 4 Fig. 4. FTIR spectra of CaHAp and CuHAp before and after sorption experiment in 0.02 M O-

Fig. 5 .

 5 Fig. 5. Mass loss, DTA, CO 2 and NO x evolvement profile curves of CaHAp, mixture of

Fig. 6 .

 6 Fig. 6. DTA and mass loss curves of CuHAp and the samples after sorption experiment in 0.02

Fig. 7 .

 7 Fig. 7. The peak area of NO x evolvement from samples after L-Ser and O-Ph-L-Ser sorption

O

  Scheme 2. P-Ser dissociation scheme[START_REF] Śmiechowski | Theoretical pKa prediction of O-phosphoserine in aqueous solution[END_REF].

Figure S3 .

 S3 Figure S3. UV absorption intensity at 200-210 nm a) for 0.01 M L-Ser solution with different Ca

Figure S7 .

 S7 Figure S7. IR spectra of : 1 -CaHAp; 2 -CaHAp after treatment in water at pH 8, 3 -CaHAp

Table 1 .

 1 Chemical composition and specific surface area of the apatites.

	Apatite	(Ca+Cu)/P,	Ca,	P,	Cu,	Specific surface area,
		Mole ratio	%	%	%	m 2 /g
	CaHAp	1.62	38.25	18.28	0	126
	CuHAp	1.62	36.10	17.54	0.91	154

  3.2.2. Release of Ca2+ and Cu2+ .The results concerning the release of Ca 2+ and Cu 2+ in aqueous solution when treating CaHAp and CuHAp by L-Ser or O-Ph-L-Ser are presented in Fig.S2. Ca 2+ and Cu 2+ release was 35 and 27 rel.%, accordingly, in L-Ser solution at pH 2. In L-Ser solution at pH 6 or 8, the Ca release is very low, especially, at pH 8, decreasing from 0.3-0.4% at pH 6 to 0.04-0.08% at pH 8. The release does not depend on L-Ser concentration. In the presence of AA the apatite solubility is a little lower than in pure water at pH˃6 (0.12%) that is in agreement with the ability of AA to give complexes with calcium ions at the apatite surface[START_REF] Spanos | Model Studies on the Interaction of Amino Acids with Biominerals: The Effect of L-Serine at the Hydroxyapatite-Water Interface[END_REF]. In contrast with Ca 2+ ions the release of Cu 2+ ions is much higher, achieving 29.3% at pH 8. As a notable difference with Ca, the Cu release is independent on the initial pH of the solution.

Substantially higher Ca release occurs in O-Ph-L-Ser solutions. As expected, due to apatite dissolution, Ca content in solutions is the highest with initial pH 2. Ca relative solubility in water with HNO 3 was 10.61 and 12.74 rel.% for CaHAp and CuHAp, respectively. It increases from 12.8 to 63.1 rel.% for both apatites with the increase of AA concentration in solution from 0.01 to 0.05 M. The Ca solubility drops to 1-3 rel.% at pH 6 and 0-5 rel.% at pH 8. The Cu release has a minimum at pH 6 and the highest values at pH 2 and 8 that is more obvious for concentrations ≥ 0.02 M. Ca release is much higher than Cu release at pH 2, but above pH 6 the Cu release exceeds Ca release. Cu solubility in O-Ph-L-Ser solutions is a little lower than in L-Ser solutions.

Table S1 .

 S1 Unit cell parameters and particles size of apatite samples before and after sorption experiments in 0.02 M O-Ph-L-Ser and L-Ser solution.

	pH 2	62(15)	9.392	6.852
	pH 6	65(15)	9.396	6.859
	pH 8	68(16)	9.390	6.854
	Sample L-Ser	Average apparent	Unit cell parameters, ±0.001Å
		size,
			a	c
	pH 6	nm 70 (19)	9.420	6.871
	CaHAp pH 8	110 (21) 76 (20)	9.413 9.431	6.860 6.878
	O-Ph-L-Ser
	pH 2	104 (23)	9.398	6.856
	pH 6	107 (21)	9.410	6.859
	pH 8	103 (20)	9.406	6.860
	L-Ser
	pH 2	108 (27)	9.444	6.885
	pH 6	98 (22)	9.440	6.882
	pH 8	104(25)	9.437	6.879
	CuHAp	68 (17)	9.394	6.855
	O-Ph-L-Ser

Acknowledgments

The current research was supported by the Estonian Ministry of Education and Research (Institutional research funding IUT33-19), CNRS, UPMC and TTU.

The authors thank Dr. Jamal Moussa for helpful discussions.


François Boulier

Marc Lefranc

François Lemaire
email: francois.lemaire@lifl.fr

Pierre-Emmanuel Morant
email: pierre-emmanuel.morant@phlam.univ-lille1.fr

Model Reduction of Chemical Reaction Systems using Elimination

There exist different schemes of model reduction for parametric ordinary differential systems arising from chemical reaction systems. In this paper, we focus on some schemes which rely on quasi-steady states approximations. We show that these schemes can be formulated by means of differential and algebraic elimination. Our formulation is simpler than the classical ones. It permitted us to obtain an approximation of the basic enzymatic reaction system which is different from those of Henri-Michaëlis-Menten and Briggs-Haldane.

Introduction

Chemical reaction systems provide a formalism to describe dynamical systems in contexts more general than the traditional chemically reacting mixtures. The dynamics of these systems is usually studied by translating them to systems of parametric ordinary differential equations by means of the mass action law. Observe that the strict setting of chemical reactions is usually not sufficient to describe many behaviours and needs to be generalized [START_REF] Horn | General mass action kinetics[END_REF]. For instance, in the process of modeling gene regulatory networks, it is common to allow unbalanced reactions, or reactions which do not conserve the mass of the reactants to describe, among other phenomenons, mRNA translation or protein degradation. See e.g. [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF].

The parametric ODE systems derived from these generalized chemical reaction systems are often too complicated for further analysis and need to be reduced. They are often overparameterized, which makes fitting methods to determine the parameters values heavy to carry out and unreliable. Moreover, for the important purpose of understanding which parts of the system contribute the most to some property of interest [START_REF] Kell | The Role of Modeling in Systems Biology[END_REF], it is preferable to deal with smaller systems.

For these reasons, model reduction of chemical reaction systems have become a central problem in their study, as explained in the survey [START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF]. Among all the existing model reduction methods (lumping, sensitivity analysis, . . . ) this paper is only concerned by the quasi-steady state approximation, which relies on the assumption that some of the chemical reactions are much faster than the other ones.

The idea of quasi-steady state approximation is simple: study the dynamics of the slow reactions, assuming that the fast ones are at quasi-equilibrium, thereby removing from the ODE system, the differential equations which describe the evolution of the variables at quasi-equilibrium. Many authors [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Okino | Simplification of Mathematical Models of Chemical Reaction Systems[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] state that carrying out rigorously this approximation is far from straightforward. We would rather say that there are different ways to perform this approximation and that it is the problem of ascertaining the domain of validity of each kind of approximation which is not straightforward.

A basic quasi-steady state approximation may very well apply if one is only interested by a basic property of the model under study (e.g. the presence of a Poincaré-Andrononv-Hopf bifurcation) possibly under some extra assumptions [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF][START_REF] Boulier | Differential Elimination and Biological Modelling[END_REF][START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF]. A more subtle one may be needed if one is interested by the timescales over which the system equilibrates or by the period and the amplitude of the oscillations, in the case of an oscillating system [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF].

This paper focuses on the quasi-steady state approximations considered in [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF], which are equivalent and which assume that reactions can be separated in only two categories: the "fast" ones and the "slow" ones (there are no, say, "very fast" or "medium" reactions). It shows how differential elimination methods [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF][START_REF] Hubert | Factorization free decomposition algorithms in differential algebra[END_REF] can be used to perform the model reduction. The presented method is brute force. A more subtle use of non differential elimination methods, based on [START_REF] Lemaire | The RegularChains library in MAPLE 10[END_REF], was implemented by the third author and will be described in a future paper. Observe however that this paper is not concerned by the problem of checking some validity conditions such as conditions C3 and C4 of [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF].

We stress the fact that using automatic elimination methods makes the task of the modeller much easier. In particular, [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF] do not simplify their reduced dynamical systems by taking into account the conservation laws. Simplifying a dynamical system modulo algebraic equations is not easy to perform by hand while it is very important since it sometimes permits to decrease the number of variables of the resulting dynamical system.

The method over the basic enzymatic reaction system

The following classical system of chemical reactions ( 1) and ( 2) describes the transformation of a substrate S into a product P under the action of the enzyme E. An intermediate complex C is produced.

E + S k1 --⇀ ↽-- k-1

C

(1)

C k2 -→ E + P (2)
There is a straightforward way to transform this chemical system into a system of four parametric ordinary differential equations. The so obtained model ( 3) is often called the initial model of the chemical system:

Ė = -k 1 E S + k -1 C + k 2 C (3a) Ṡ = -k 1 E S + k -1 C (3b) Ċ = k 1 E S -k -1 C -k 2 C (3c) Ṗ = k 2 C (3d)
The hypothesis which permits to perform all the model reductions is that the reversible reaction [START_REF] Bellon | Biochimie structurale : introduction à la cinétique enzymatique[END_REF] is much faster than the reaction [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF]. With other words, one assumes that

k 1 , k -1 ≫ k 2 .
Two different reductions of this chemical system were given in the early twentieth century by Henri, Michaëlis and Menten [START_REF] Henri | Lois générales de l'Action des Diastases[END_REF][START_REF] Michaelis | Die kinetik der invertinwirkung[END_REF] on the one hand, Briggs and Haldane [START_REF] Briggs | A note on the kinetics of enzyme action[END_REF] on the other hand.

The Henri-Michaëlis-Menten and Briggs-Haldane reductions

In all cases, the reduced model, which describes the evolution of the product concentration from the substrate concentration, writes:

Ṗ (t) = V m S(t) K + S(t) , Ṡ(t) = - V m S(t) K + S(t) (4)
but with different assumptions and different values for the parameter K.

Both reductions rely on a few extra assumptions i.e. (the 0 subscript indicating initial concentrations) S 0 ≫ E 0 , P ≃ 0, C 0 = 0 [1, page 9]. In the case of the Henri-Michaëlis-Menten reduction, one assumes equation ( 1) is constantly at equilibrium (the pre-equilibrium condition), which means k 1 C ≃ k -1 E S. Thus Ṡ is assumed to be small. Performing some extra computations, one is led to V m = k 2 E 0 and K = k-1 k1 • In the case of the Briggs-Haldane reduction, one assumes that the intermediate complex does not vary much i.e. that Ċ ≃ 0 which implies

k 1 E S ≃ (k -1 + k 2 )
C (the quasi-stationary state hypothesis). Doing some further computations, one is led to V m = k 2 E 0 (as in the Henri-Michaëlis-Menten case) and K = k-1+k2 k1 •

Reduction using differential elimination

The first idea that a casual reader would have for performing such reductions could consist in translating the "≃" into a "=". In the Briggs-Haldane case, the hypothesis Ċ ≃ 0 would be treated by enlarging the initial model equations with the new equation Ċ = 0. However, this equation implies that C is a constant, which implies that S is a also a constant. The reduced model so obtained can certainly be considered as useless.

The right way to proceed consists in separating the "fast" variables from the "slow" ones and in transforming the ODE defining the evolution of each fast variable ẋ = rhs by the algebraic equation 0 = rhs. But what is a "fast" variable ? An intuitive definition would be: a variable is fast if it is involved in a fast chemical reaction. Over some examples [START_REF] Vilar | Mechanisms of noise-resistance in genetic oscillators[END_REF] this definition is sufficient. However, it would clearly lead to a useless reduced model over our example since E, S and C would then be considered as fast variables and one would end up with the reduced model Ṗ = Ṡ = 0. As explained in [2, page 3502], a variable x such that ẋ depends on a mixture of slow and fast reactions should not be considered as slow. The fast and slow variables of a system are not necessarily the model variables but functions of them. They can be computed by a rigorous two timescales analysis, as shown by [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF] but, as stated by [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF], this is not always necessary.

The reduction that we present here is equivalent to those of [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] but differs in its presentation. As for the Henri-Michaëlis-Menten reduction, one uses the pre-equilibrium approximation on the first reaction but this is (for the moment) the only assumption we make:

k 1 E S = k -1 C (5)
As for the initial model, one translates the two chemical reactions as parametric ordinary differential equations expressing the evolution of the concentrations of E, S, C and P . However one handles the reactions which are assumed to be at equilibrium (the fast reactions) in a different manner than the other reactions (the slow reactions): the variation of the concentration of each chemical species is the sum of the contributions of each reaction. The contribution of each slow reaction is derived from the law of mass action as usual. The contribution of each fast reaction is represented by a new variable.

The underlying idea is that the dynamic of the fast reaction is (for the moment) considered as unknown. This means that one adds a degree of freedom for each fast reaction. This freedom in fact allows the system to stay on the equilibrium conditions. Performing elimination on the extra unwanted variables, one gets a set of differential equations in the reactants only. Over the example, the contribution of the fast reaction is denoted F 1 (for fast reaction 1) and the contribution of the slow reaction is k 2 C:

Ė = -F 1 + k 2 C (6a) Ṡ = -F 1 (6b) Ċ = F 1 -k 2 C (6c) Ṗ = k 2 C (6d)
Using elimination on ( 5) and ( 6)

with the ranking R = [F 1 ] ≫ [P, C, E, S] [16, chap. I, §8], one eliminates the unknown F 1 .
Elimination is performed below using the MAPLE diffalg package. The parameters are put in the base field F of the equations to avoid discussions over their values (they are then considered as algebraically independent). An inequation is added to avoid considering the degenerate case of C(t) being the zero function. The output is pretty printed. with(diffalg): rules := rewrite_rules( Ids [START_REF] Bellon | Biochimie structurale : introduction à la cinétique enzymatique[END_REF] );

sys := [ E[t] -(-F1 + k2*C), S[t] -(-F1), C[t] -(F1 -k2*C), P[t] -k2*C, k1*E*S -
F 1 = k 2 k 1 E S (k 1 S + k -1 ) k -1 (k 1 S + k 1 E + k -1 ) , Ė = k 1 2 E 2 k 2 S k -1 (k 1 S + k 1 E + k -1 ) , Ṗ = k 2 k 1 E S k -1 , Ṡ = - k 2 k 1 E S (k 1 S + k -1 ) k -1 (k 1 S + k 1 E + k -1 ) , C = k 1 E S k -1
In the above system, one gets the value of F 1 . The reduced model is made of three parametric ordinary differential equations. It is exactly that obtained by [25, page 2327, (33), denoting A, B, C, D, κ for S, E, C, P, k 1 /k -1 ].

Refinements

One can simplify further the reduced model by replacing the two parameters k 1 and k -1 by a single parameter K = k -1 /k 1 . Indeed, k 1 and k -1 only appear in the equation ( 5) for which one can substitute K = k -1 /k 1 . This simplification makes sense: the values of k -1 and k 1 describe the speed of the fast reaction. Since this reaction is assumed to be constantly at the equilibrium, only the ratio of the two constants matters. Continuing the example, one gets: map (normal, subs (km1=K*k1, rules));

" F1 = k2 E S (S + K) K (S + E + K) , Ė = k2 E 2 S K (S + E + K) , Ṗ = k2 E S K , Ṡ = - k2 E S (S + K) K (S + E + K) , C = E S K # .
The second simplification consists in using the two conservations laws below (that one can automatically deduce from the system (6) by means of linear algebra):

E + C = E 0 + C 0 (7a) S + C + P = S 0 + C 0 + P 0 (7b)
plus the hypothesis C 0 = P 0 = 0 (as for the classical reductions). One can then eliminate E from the ODE describing the evolution of S. This can be done by means of elimination. For instance, one can rerun the above commands, enlarging the field F with the two new symbols E 0 and S 0 and the list sys with the two equations E + C = E 0 and S + C + P = S 0 . Performing the same elimination as above and denoting k 2 E 0 = V m as usual, one eventually gets:

Ṡ = - V m S (K + S) K E 0 + (K + S) 2 • (8)
As far as we know, this formula is new. If one neglects the term K E 0 in equation ( 8), one gets equation [START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF]. Thus equation ( 8) and ( 4) are almost equal when the term K E 0 is negligible, for example when S ≫ E 0 . What happens when S and E 0 are in the same range, a phenomenon likely to occur in vivo according to [10, §2.1.1] ? Equation ( 8) turns out to be more precise, Figure 1 shows.

The algorithm

The main algorithm DifferentialModelReduction is composed of two main steps: the transformation of the list of reactions into a dynamical system, and the elimination part. Its main interest is that it is fully automatic and performs automatic case splitting when needed.

The DifferentialModelReduction algorithm described in section 3.2 relies on a few auxiliary functions described in section 3.1. In all algorithms presented in this section, the parameter L denotes a list of chemical reactions, and X denotes a vector of all chemical species of L (i.e. all the reactants and products involved in L). Except in the function DifferentialModelReduction, the Figure 1. When the initial enzyme concentration is close to that of the substrate, the curve S(t) computed from the initial model (3) (red, coming from the top of the diagram) is better approximated by formula (8) (green) than by formula (4) (blue). The two reduced formulas were integrated for an initial value picked on the non reduced curve, after the transient step, at t = 0.005. Parameters values were borrowed from [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF]. order of the elements of X is not important. As explained later in section 3.2, the order of X encodes the chemical species to be eliminated.

Each reaction is represented by a list of reactants, a list of products, a boolean flag indicating if the reaction is considered fast, a boolean flag indicating if the reaction is reversible, and one (or two) rate constant(s) (depending on whether the reaction is reversible or not). Reactants and products are represented by lists with two elements: a symbol and a stoechiometric coefficient. For example, the two reactions (1) and (2) in section 2 could be coded in MAPLE, in the humanly readable format:

[ reactants=[ [E,1], [S,1] ], products=[ [C,1] ], rateConstants=[k1, km1], fast=true, reversible=true ]; [ reactants=[ [C,1] ], products=[ [E,1], [P,1] ], rateConstants=[k2], fast=false, reversible=false ]
Observe that, in our algorithm, a reversible reaction is not equivalent to two irreversible reactions. The task of indicating that a reaction is reversible is left to the user. It could be made automatic by computing the rank of a stoechiometry matrix as [25, page 2323].

A basic way to define slow and fast reactions consists in splitting the set of rate constants into two sets: the large ones and the small ones. A reaction is considered as fast if it involves at least one large rate constant. Stricly speaking, the notion of fast reaction should take into account the concentrations.

The auxiliary functions

In this section, one denotes a reaction as

a j X j k1 -→ b j X j if R is irreversible, otherwise a j X j k1 -⇀ ↽ - k2 b j X j .
This formulation involves all the chemical species present in the list of reactions, by setting a j = 0 (resp b j =0) when X j is not a reactant (resp a product). The function Equilibria. It returns a list of algebraic relations. It simply converts each fast reaction into a single algebraic equation coding the pre-equilibrium. On our example in section 2, it simply returns [k 1 ES -k -1 C], which is the equation ( 5).

Algorithm 1 Equilibria(L, X)

Input: a list of reactions L, a vector X of all reactants and products involved in L Output: the list of equilibrium relations for the fast reactions 1: equis := {} 2: for each fast reaction R of L do

3: denote R as a j X j k1 -→ b j X j if R is irreversible, otherwise a j X j k1 -⇀ ↽ - k2 b j X j 4:
if R is reversible then equis := equis ∪ {eq} 10: end for 11: return equis The function ListReactionsToODE. It returns a dynamical system Ẋ = H(X, Θ, F ), where Θ denotes the set of all rate constants in the list of reactions, and F denotes the set of the extra variables F i introduced for the fast reactions. On our example in section 2, if L is the list composed Algorithm 2 ListReactionsToODE(L, X)

5: eq := k 1 X j aj -k 2 X j bj 6:
Input: L = [R 1 , . . . , R p ]
is a list of reactions, X = t (X 1 , . . . , X n ) is a vector of all reactants and products involved in L Output: a dynamical system Ẋ = H(X, Θ, F ) describing L 1: {let us first construct the vector of rate V } 2: let V = t (v 1 , . . . , v n ) 3: for each reaction l i do 4:

denote l i as a j X j k1 -→ b j X j if R is irreversible, otherwise a j X j k1 -⇀ ↽ - k2 b j X j 5:
if l i is a fast reaction then 6:

v i := F i 7: else 8:
if l i is irreversible then 9:

v i := k 1 X j aj 10:
else 11:

v i := k 1 X j aj -k 2 X j bj 12:
end if m ij = b j -a j for each 1 ≤ j ≤ n 20: end for 21: return Ẋ = t M.V of the two reactions (1) and (2) (in that order), and if X = t (E, S, C, P ), then the vector V and the matrix M computed inside ListReactionsToODE are

V = F 1 k 2 C M = -1 -1 1 0 1 0 -1 1 Thus ListReactionsToODE returns Ẋ =     Ė Ṡ Ċ Ṗ     =     -1 1 -1 0 1 -1 0 1     F 1 k 2 C =     -F 1 + k 2 C -F 1 F 1 -k 2 C k 2 C    
which is exactly the system (6). The ConservationLaws function. It aims at finding linear combinations of the X i which are constant along solutions. To compute them, one just needs to compute a basis of the kernel of M (i.e. a basis of the solutions of M.Y = 0). This is a linear algebra problem over Q. The argument is straightforward: if y satisfies M.y = 0, then t y. Ẋ = t (M.y).V = 0, so y.X is constant on any solution, which implies y.X is a conservation law.

Algorithm 3 ConservationLaws(L, X)

Input: a list of reactions L, a vector X of all reactants and products involved in L Output: a list of (linear) conservation laws 1: denote M the stoechiometric matrix as defined in ListReactionsToODE 2: compute a basis [y 1 , . . . , y s ] of M.Y = 0 3: return the list [y 1 .X, . . . , y s .X]

On the example of section 2, a basis of the kernel of M is for example t (-1, 1, 0, 1) and t (1, 0, 1, 0). This leads to -E + S + P = -E 0 + S 0 + P 0 and E + C = E 0 + C 0 . Although different to the two conservation laws [START_REF] Boulier | Computing canonical representatives of regular differential ideals[END_REF], one easily checks that they are equivalent. If instead, one takes the other basis t (0, 1, 1, 1) and t (1, 0, 1, 0), one exactly obtains the conservation laws [START_REF] Boulier | Computing canonical representatives of regular differential ideals[END_REF].

The DifferentialModelReduction algorithm

The function DifferentialModelReduction. It transforms a list of reactions into a list of dynamical systems encoding the reduced models. The chemical species to be eliminated are to be put first in the vector X. The number of chemical species that can be eliminated depends on several criteria: the more fast reactions and conservation laws there are, the more chemical species get eliminated. However, this number also depends on the form of the equations themselves, as shown in section 4.1.

After building the system S of equations, DifferentialModelReduction builds the ranking R to be used by Rosenfeld-Gröbner over the system S. The ranking [F ] ≫ [X] means that all derivatives of F are greater than those of X, the derivatives of F are ordered w.r.t. an orderly ranking, and the derivatives of X also.

The algorithm Rosenfeld-Gröbner performs the so called differential elimination. It returns a list of systems of polynomial differential equations called characteristic sets. Each characteristic set C i defines some differential ideal A i in the differential polynomial ring R which contains the model equations. Using the set C i , one is able to perform computations (such as the normal form computation) modulo the ideal A i . The italicized words above come from the differential algebra [START_REF] Ritt | Differential Algebra[END_REF][START_REF] Kolchin | Differential Algebra and Algebraic Groups[END_REF].

Applied to any differential polynomial p and C i , the NormalForm algorithm returns a rational fraction which is a canonical representative of the residue class of p in the factor ring R/A. See [START_REF] Boulier | Computing canonical representatives of regular differential ideals[END_REF] or [8, Figure 1] for a presentation of the NormalForm algorithm 1 . In our case, the goal of NormalForm is straightforward: it consists in expressing each derivative Ẋi or each unknown F j Algorithm 4 DifferentialModelReduction(L, X) Input: a list of reactions L, a vector X = t (X 1 , . . . , X n ) of all reactants and products involved in L Output: a list of reduced dynamical systems in the variables X 1: Sys := ListReactionsToODE(L, X) 2: Cons := ConservationLaws(L, X) 3: Equi := Equilibria(L, X) 4: S := Sys ∪ Cons ∪ Equi 5: R the ranking [F ] ≫ [X] where F is the set of the unknowns F i introduced for each fast reaction 6: [C 1 , . . . , C t ] := Rosenfeld-Gröbner(S, R) 7: if NormalForm(F i , C k ) is not a rational fraction in the variables X (for all F i and all 1 ≤ k ≤ t) then FAIL.

8: return the list [S 1 , . . . , S t ]
where each

S i is equal to [ Ẋ1 = NormalForm( Ẋ1 , C i ), . . . , Ẋn = NormalForm( Ẋn , C i )]
as a rational fraction in the variables X only. If this goal is reachable, the normal form function performs it, thanks to the rankings properties.

Rosenfeld-Gröbner returns a list of characteristic sets because it is a case splitting algorithm, so it may need to consider different cases (say) X j = 0 and X j = 0 for some j. Such splittings might introduce some spurious outputs, since the case X j = 0 is usually not very interesting. Indeed in section 2, the case C = P = E = S = 0 is not interesting. However, those splittings are needed when fast irreversible reactions are involved. For example, if one assumes that the reaction A + B k -→ C is fast, then the pre-equilibrium hypothesis yields kAB = 0, which implies A = 0 or B = 0. Those two cases need to be treated differently since the dynamic in each case is different.

Over our example in section 2, the function DifferentialModelReduction (up to the renaming K = k -1 /k 1 ) performs the same computations as in sections 3.2 and 2.3.

Failures

Strictly speaking, it may happen that some unknowns F i cannot be expressed as rational fractions of the variables X. In this case, our algorithm fails. Examples for which our method fails include many examples which make no sense (e.g. fast unbalanced reactions such as A → A + B). A meaningful example for which our method does not readily apply is given by the Robertson's example [12, chap. IV.1] when the fast and the very fast reactions are considered both fast.

Translating our technical condition in terms of properties of the chemical reaction systems is an interesting problem, left for a future paper. Observe that our condition plays a role equivalent to the "kinetic" linear independence condition combined to the singularity testing for algebraic equations described in [25, pages 2324-2326] but our method seems to be more straightforward: it can be directly formulated in the setting of algebraic geometry (ideal theory) instead of a matrix inversion problem over a residue class ring [25, eq. ( 27)].

Sufficient conditions for the two methods to apply are described in [25, page 2326] (e.g. case of fast reactions with pairwise disjoint sets of reactants and products).

Optimisations

First, as said in section 2, one can decrease the number of rate constants in the reduced model, when there are reversible fast reactions. Indeed, for each fast reversible reaction

a j X j k1 -⇀ ↽ - k2 b j X j ,
the function Equilibria could simply compute eq := X j aj -K X j bj where K is a new symbol equal to k 2 /k 1 . This leads to simpler equations. To handle this algorithmically, one would need to modify the function Equilibria so that it also returns a description of the new K symbols.

Second, although the differential elimination is conceptually simpler, it appears that the elimination process is mainly a linear algebra problem over some residue class ring. This can be achieved by using regular chains and regularity tests, as can be done with the RegularChains [START_REF] Lemaire | The RegularChains library in MAPLE 10[END_REF] package (and its subpackage MatrixTools) in MAPLE. Here is some insight of this claim (another paper will discuss that point).

Over the example of section 2, differentiate (5):

k 1 ĖS + k 1 E Ṡ = k -1 Ċ (9)
Putting ( 6) and ( 9) together, one gets a system of five linear equations in F 1 , Ė, Ṡ, Ċ and Ṗ , over the residue class of K[E, S, C, P, k 1 , k -1 , k 2 ] by [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF]. This system can be solved w.r.t. F 1 . Taking conservation laws into account preserves the linearity of the problem.

Examples

In section 4.1, we carry out a basic polymer degradation system. We checked that our method applies to more complicated systems such as a polymerisation cascade or a variant of the ozone decomposition system (dropping temperature considerations) [25, page 2328].

Polymer degradation

One considers a simple dimerisation between a protein P and itself. One also assumes that the protein is degraded. The dimerisation is assumed to be fast, the degradation being slow.

2P k1 --⇀ ↽-- k-1 P 2 (10a) P k2 -→ ∅ (10b)
Let us introduce K = k 1 /k -1 and take X = t (P, P 2 ). If L represents the two reactions (10a) and (10b), the call to DifferentialModelReduction(L, X) performs the following computations. The list Sys simply consists of

Ṗ = -2F 1 -k 2 P (11a) Ṗ2 = F 1 (11b)
The

F 1 = -2 k 2 P 2 (4 KP -1) 16 KP 2 -1 , Ṗ2 = -2 k 2 P 2 (4 KP -1) 16 KP 2 -1 , P 2 = P 2 K
Computing the normal forms of Ṗ and Ṗ2 , one gets:

Ṗ = - k 2 P 2 (4 KP -1) KP (16 KP 2 -1) (12a) Ṗ2 = -2 k 2 P 2 (4 KP -1) 16 KP 2 -1 (12b)
In this case, one could not eliminate P since it still appears in the equation (12b). However, over this simple example, one could eliminate P by using the (non algebraic) relation P = P 2/K. This would lead to: Ṗ2 = -

2k 2 P 2 4 √ KP 2 + 1 (13)
However, this cannot be generalized easily (since an algebraic equation cannot be symbolically solved when its degree is greater than 5). If one changes the order and takes X = t (P 2 , P ), one obtains a different system:

Ṗ2 = - 2k 2 KP 2 4KP + 1 (14a) Ṗ = - k 2 P 4KP + 1 (14b)
In that case, P 2 is eliminated.

Equivalence of our method with [25, 23, 2]

Our method is essentially the same as that of [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF]. The main differences are: we do not introduce any ǫ but the unknown F i instead, we perform differential elimination instead of linear algebra and we directly handle reversible reactions instead of splitting them into two irreversible reactions (this affects the rate vector in the ListReactionsToODE function).

The above comments also apply to [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF] since this paper is a basis for [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF]. The paper [START_REF] Van Breusegem | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF] considers more general systems with external supply and withdrawal of species. For this reason, they do not obtain our formula (8) which is a consequence of conservation laws which do not exist in this general setting. They provide a strong theoretical basis for all methods. They transform the initial system into the so called two-time scale standard form by expliciting the true slow variables in order to apply the Tikhonov theorem. Our method (as well as that of [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF]) is equivalent but more direct since it avoids this change of coordinates.

Our method is also essentially the same as that of [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF]. The paper [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] presents an approximation scheme using the so called prefactors. Though we do not introduce this concept, our algorithm DifferentialModelReduction gives exactly the same reduced systems as [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] for the examples listed in [2, tables 1 and 2]. A major difference is that our method is fully algorithmic: it does not need to explicit the true slow variables and relies on elimination to perfom the computations. On the opposite, [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] does need to explicit the true slow variables (which might actually be computed automatically, following [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF] for example). Computations are then performed by hand.

In order to illustrate the equivalence of our reduction method with that of [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF], let us consider a simplified version of the example provided in [2, tables 1]: the case of the simple dimerisation given in section 4.1. Mutatis mutandis, the argumentation of [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] writes as follows.

One assumes that the dimerisation is fast compared to the degradation. This means that one has P 2 = K P 2 . Let us introduce n P the total number of proteins P , which is simply n P = P +2 P 2 (since a dimer P 2 contains two proteins P ). This new variable is in fact a slow variable, because this quantity does not change even when the fast reaction occurs. The total degradation on P is simply the degradation on P , since P 2 is not degraded. Therefore one has ṅP = -k 2 P . By using the three previous equations, one can deduce Ṗ = -k2 P 4 K P +1 by doing the following computations. One has ṅP = -k 2 P . One also has n P = P + 2P 2 , which implies ṅP = Ṗ + 2 Ṗ2 . Using P 2 = KP 2 , one gets using a differentiation Ṗ2 = 2KP Ṗ . Therefore, one has -k 2 P = Ṗ + 4KP Ṗ , which yields Ṗ = -k2P 4KP +1 Indeed, the approach developed in [START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] introduces the three equations

P 2 = KP 2 (15a)
n P = P + 2 P 2 (15b) ṅP = -k 2 P (15c) whereas our approach introduces:

P 2 = KP 2 (16a) Ṗ = -2F 1 -k 2 P (16b) Ṗ2 = F 1 (16c)
Although based on different variables, systems [START_REF] Kell | The Role of Modeling in Systems Biology[END_REF] and ( 16) are in fact very similar since they both imply:

P 2 = KP 2 (17a) Ṗ + 2 Ṗ2 = -k 2 P (17b)
Indeed equation (17b) is obtained either by combining (16b) and (16c), or by combining (15b) and (15c). Since system [START_REF] Lemaire | The RegularChains library in MAPLE 10[END_REF] implies Ṗ = -k2P 4KP +1 , both systems ( 15) and ( 16) also do.

Conclusion

We have presented an algorithmic method for reducing systems of chemical reactions using elimination and we have shown the relationship between the approaches of [START_REF] Vora | Nonlinear model reduction of chemical reaction systems[END_REF][START_REF] Bennet | Transient Dynamics of Genetic Regulatory Networks[END_REF] and ours.

The main interest of our method is that it is fully automated (once the set of fast reactions has been chosen). Moreover the setting of our algorithm could easily be integrated in software based on the SBML format since recent versions of SBML (level 2 version 3 for example) include the flags reversible and fast that we need in order to apply our method. Still in this computational context, the reduction of the number of rate constants based on symmetries [START_REF] Sedoglavic | Reduction of algebraic parametric systems by rectification of their affine expanded lie symmetries[END_REF] could also be used.

The replacement of the differential elimination step by linear algebra methods over residue class rings is left for a future paper as well as a formulation of the failure conditions in chemical terms.

The application of DifferentialModelReduction to the systems of [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Morant | Oscillateurs génétiques simples[END_REF] and the analysis of the reduced models are also left for a future paper. In the case of [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF], the new reduced system is more complicated to handle than the previous one but proves to be a better approximation of the initial model it comes from.

  km1*C ]: F := field_extension (transcendental_elements=[k2,k1,km1]): R := differential_ring (ranking=[ [F1], [C,E,P,S] ], derivations=[t], field_of_constants=F): Ineqs := [ C ]: Ids := Rosenfeld_Groebner (sys,Ineqs,R); Ids := [characterizable]

  list Cons is empty, and Equi is [KP 2 -P 2 ]. Let us detail the computations using diffalg: Sys := [ P[t] -(-2*F1 -k2*P), P2[t] -F1 ]; Cons := []; Equi := [ K*P^2-P2]; F := field_extension (transcendental_elements=[K,k2]): R := differential_ring (ranking=[ [F1], [P,P2] ], derivations=[t], field_of_constants=F): Ineqs := [ P ]: Ids := Rosenfeld_Groebner ( [op(Sys),op(Cons),op(Equi)] ,Ineqs,R); rules := rewrite_rules( Ids[1] );

In very rare cases, the NormalForm algorithm may split the characteristic set C i as finitely many smaller characteristic sets, refining thereby the list returned by Rosenfeld-Gröbner.