Sleiman Rita
email: rita.sleiman@ensait.fr

Tran Kim-Phuc

Thomassey Sébastien

Long Term Demand Forecasting System for Demand Driven Manufacturing

Keywords: Demand Driven Manufacturing, Sales forecasting, Demand forecasting, Historical sales correction

Demand-Driven Manufacturing (DDM) is the solution that most companies are heading to in our days. Although this strategy consists of producing goods based on what consumers demand, companies should also rely on accurate forecasting systems to prepare their production chain for such an operation by supplying enough raw material, increasing production capacity to fit the desired demand, etc.… However, due to the fact that most companies have been relying on massive production, most sales forecasting systems usually used rely on sales data of previous years that, not only contain the actual demand, but takes into consideration the marketing strategy effects like massive promotions. Hence, the resulting forecasts do not mainly reflect consumers' demand. For this reason, a switch to demand forecasting, instead of sales forecasting, is essential to ensure a good transition to DDM. This paper proposes an artificial intelligence based demand forecasting system that aims to determine "potential sales", mainly reflecting consumers' demand, by correcting historical sales data from external variables' effects. A comparison with other sales forecasting models is performed and validated on real data of a French fashion retailer. Results show that the proposed system is both robust and accurate, and it outperforms all the other models in terms of forecasting errors.

Introduction

For many years, the market of consumer goods relies on overconsumption with mass production and aggressive marketing campaigns. This approach requires the companies to focus on low cost production with off-shore sourcing and involves both high volume purchases and long lead time [START_REF] Koren | Mass Production and Lean Manufacturing Book: The Global Manufacturing Revolution: Product-Process-Business Integration and Reconfigurable Systems[END_REF].

Due to the needs for more personalized products and environmental awareness of the consumers, this strategy reaches its limits in many industrial and commercial sectors [START_REF] Hu | Evolving paradigms of manufacturing: from mass production to mass customization and personalization[END_REF]. Consequently, many companies move toward a Demand Driven Manufacturing (DDM) strategy [START_REF] Lebovitz | The journey toward demand driven manufacturing[END_REF]. DDM strategy aims to produce goods based on consumers' demand instead of sales or target forecasts [START_REF] Tiedemann | Demand-driven supply chain operations management strategies -a literature review and conceptual model[END_REF]. According to Gartner, nearly 90% of production industries which are not demand-driven, are planning to be [START_REF]What is Demand-Driven Manufacturing[END_REF]. This is mainly because of the benefits a DDM approach could offer, such as increasing profitability, reducing wastefulness, and most importantly improving customers' satisfaction.

Even though a DDM strategy relies on flexible and reactive production and supply chain, companies still require accurate forecasting systems to anticipate their raw material supplies, production capacity, … and meet the consumer expectations in terms of lead time. However, current sales forecasting techniques do not reflect exactly consumer's demand since they rely on historical data including the constant discounts or commercial offers. Moreover, in a DDM approach, only goods that are going to be sold are produced, which implies the application of a totally different marketing strategy. Hence, in order to ensure a successful transition to the DDM strategy, companies should switch from a sales forecasting system to a demand forecasting system.

Recently, Artificial Intelligence AI techniques have presented great improvements in many fields and especially in forecasting [START_REF] Raza | A review on artificial intelligence based load demand forecasting techniques for smart grid and buildings[END_REF][START_REF] Thomassey | Introduction: Artificial Intelligence for Fashion Industry in the Big Data Era[END_REF]. In fact, with the AI techniques and the heterogeneous data sources available nowadays, new forecasting systems can be developed to determine consumer demand by modeling the impact of the different exogenous variables on sales.

This paper proposes a demand forecasting system that aims to correct historical sales, using AI techniques, by removing the influences caused by external factors, such as promotions, before moving to the prediction process. In fact, by modeling the impact of the previous discount strategy, historical sales could be corrected, presenting then potential sales that mainly reflect consumers' demand. This results in an accurate and robust forecasting model which can be an additional tool to help companies better react to the rapidly changing market, and meet the demand of materials, warehouses, etc… for the production process.

The following parts will be structured as follows. A brief literature review about existing sales forecasting systems is presented in Sect. 2. Section 3 describes the proposed methodology in detail. Results and discussion are presented in Sect. 4. The conclusion is presented in the last section.

Literature Review

Retail companies, based on a mass production strategy, generally rely also on sales forecasting systems in order to optimize their production and supply chain. Hence, sales forecasting has been always studied by researchers in the literature, which presents a wide variety of sales forecasting systems, adapted to different contexts [START_REF] Choi | Fast fashion sales forecasting with limited data and time[END_REF][START_REF] Chang | A hybrid model by clustering and evolving fuzzy rules for sales decision supports in printed circuit board industry[END_REF][START_REF] Chen | Sales forecasting by combining clustering and machine-learning techniques for computer retailing[END_REF].

The most commonly used methods for sales forecasting are statistical time series analysis tools that primarily seek to identify trends and seasonality. Among the existing statistical methods, "Auto-Regressive Integrated Moving Average" ARIMA and "Seasonal Auto-Regressive Integrated Moving Average" SARIMA are widely used for sales forecasting [START_REF] Liu | Sales forecasting for fashion retailing service industry: a review[END_REF]. Similarly, SARIMAX is an extension of SARIMA when explanatory variables are available. Moreover, deep learning has recently shown promising results in nonlinear sequence learning problems. Notably, Recurrent Neural Network (RNN) and Long Short Term Memory (LSTM) networks are popular deep learning techniques and have outperformed popular machine learning methods for time series forecasting [START_REF] Fischer | Deep learning with long short-term memory networks for financial market predictions[END_REF]. RNN and LSTM networks, unlike other neural networks, have the property of retaining information over time.

The methods already mentioned are likely to be adequate to learn seasonality in a time series. However, in the context of sales forecasting, there are many exogenous variables that can influence the predictions. Thus, several hybrid methods that can take advantage of both time series analysis and explanatory variables are proposed in the literature [START_REF] Liu | Sales forecasting for fashion retailing service industry: a review[END_REF]. Punia et al. [START_REF] Punia | Deep learning with long short-term memory networks and random forests for demand forecasting in multi-channel retail[END_REF] proposed a new forecasting method that combines LSTM, to model the time aspect of the input series, and Random Forest (RF), to predict the impact of explanatory variables not detected by LSTM. The obtained results showed that their method outperformed other methods, ARIMAX, LSTM and RF. Another more interpretable working methodology is to correct the sales history by removing the impact of explanatory variables before making the prediction of the time series, as proposed by Thomassey et al. [START_REF] Thomassey | A short and mean-term automatic forecasting system--application to textile logistics[END_REF].

Methodology

The main goal of this paper is to build a demand forecasting based DDM approach as illustrated in Fig. 1. Using the demand forecasts, production managers can better manage the supply chain by anticipating the required raw materials and adapting their production capacity to meet actual consumers' demand.

The approach adopted in this paper was inspired by the forecasting system proposed in [START_REF] Thomassey | A short and mean-term automatic forecasting system--application to textile logistics[END_REF] in which a fuzzy logic system is optimized to predict the influences of the explanatory variables on sales. The work done in this paper uses a similar approach for correcting historical sales (year ∈ [1, . . . , n]) by removing the influence of external variables before predicting sales for year n + 1. The following parts describes in details the proposed methodology as well as the learning process.

Global Forecasting System

The proposed forecasting system, presented in Fig. 2, is composed of three main stages:

• The influences of explanatory variables are removed from historical sales data using the Influence Capturing System ICS, • The resulting historical potential sales are used to make predictions for next year's potential sales, • The final sales predictions are obtained by adding the influences of explanatory variables corresponding to the next year to the previous predictions using ICS -1 .

Learning Process

The Influence Capturing System ICS is designed to predict, using the explanatory variables, a coefficient that has a role of correcting the historical sales, generating then a series of historical potential sales. We used a simple neural network NN (a perceptron) to predict these correction coefficients, based on the assumption that a NN can learn to generate efficient and sophisticated rules for prediction.

Thereafter, for predicting the potential sales of the next year, we applied the seasonality average, which is one of the simplest time series forecasting methods. We selected this basic model in order to check the ability of the model in the optimization of the correction coefficients, without having additional parameters to adjust, which increases the complexity in the overall learning process.

Fig. 3. Learning process

Figure 3 above presents the learning process of the proposed system in detail. In fact, the NN must be optimized in order to predict the most reliable and efficient correction coefficients. To do so, we used a genetic algorithm GA where each chromosome is formed by real genes representing the weights of the neural network. The learning loop is composed of the following important steps:

• Using the explanatory variables as inputs, the NN predicts the corresponding correction coefficient (cc), which is considered to be a multiplicative coefficient that quantifies the impact of exogenous variables on sales, • The Influence Extraction function IE calculates the potential sales as input to the prediction model, • The V-curve is the predicted life curve obtained from the potential sales (without influences of the explanatory variables), • The influences of the explanatory variables are re-added to the obtained V-curve, using IE -1 , to perform the final prediction, • The GA's fitness function is calculated from the sales thus predicted and the actual sales by calculating the root mean square error RMSE:

fitness = 1/RMSE (1)
• The GA sets new weights for the neural network to increase the value of the fitness function, • The optimal parameters of the neural network that allows the determination of the optimal life curve V* (corrected historical sales) are obtained when the maximum number of iteration of the GA is reached.

Once the optimal ICS model is obtained, the final forecast is calculated by adding the influence of the explanatory variables of the future season (the 52 weeks to be predicted), to the optimal life curve V*.

Results and Discussion

Data Used

The proposed forecasting system was applied to real world data of a French fashion retailer. It was tested to perform a one year sales prediction (52 weeks) of a T-shirt's family of products. We used three years of historical sales data (2016, 2017, 2018) for the learning process, and one year, 2019, for testing.

In order to test the ability of the proposed system to capture the influence of the explanatory variables on the sales, those variables have to be wisely selected and must be the most significant ones. Therefore, since obtaining the influence of promotions on sales was important to our industry partners, the explanatory variables chosen were:

• The average percentage of discount per family of products, • The type of commercial operation OP applied (no OP, winter/summer sales of long duration, or 10-days winter/summer sales of short duration), • The season (fall, winter, spring, summer),

• The remaining duration of the sale period (100% at the beginning of the sale and 0% at the end) to locate at which phase of the sales period we are.

Prediction Results

Figure 4 shows, for the year 2019, the optimal V-curve representing the prediction without the influence of exogenous variables, the predicted sales after reinjection of these influences, the corresponding predicted corrective coefficients for each week, as well as a visualization of some of the explanatory variables used. The results obtained demonstrate the effectiveness of the proposed method as well as the ability of the neural network, combined with GA, to model the influences of the explanatory variables. It is clear from Fig. 4 that during winter/summer sales (shown in red dashed circles), the influence of the explanatory variables is more visible than in normal weeks. The correction coefficient is at its highest at the beginning of the winter/summer sales, and decreases continuously throughout the promotions period: we can clearly see in the predictions for the year 2019 that the peak in weeks 26 and 27 is caused by the beginning of the summer sales (corrective coefficient, in this case, is equal to 90%). This correction is of less importance during 10-days sales (shown in black dashed circles), which is logical since during this type of operation, promotions are not applied on the entire collection which explains the lower peaks in those periods. We evaluated our proposed methodology by comparing it to three other forecasting methods in terms of Root Mean Squared Error (RMSE), Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE). We used for comparison two of the classical time series forecasting methods, SARIMA and SARIMAX. In addition, we 1) show that the proposed methodology yielded better results when applied to 2019 test data, by having the lower forecasting errors for the three of the metrics used. Figure 5 shows a comparison of the prediction results obtained by the different methods mentioned above.

Conclusion

This paper presents a demand forecasting system that is able to correct historical sales by eliminating the effect of marketing strategies such as promotions. The proposed system plays a vital role for companies when shifting to DDM strategy. It helps them optimizing the supply chain by adapting their resources to the demand. It was shown that it is able to perform well compared to other existing sales forecasting systems due to the fact that it relies on potential sales to perform the predictions. This system can be improved with the help of experts in fields like marketing, fashion, … in order to add additional features to the learning process for obtaining potential sales that better reflect consumers' demand.

Fig. 1 .

 1 Fig. 1. Demand forecasting based DDM approach

Fig. 2 .

 2 Fig. 2. Global forecasting system

Fig. 4 .

 4 Fig. 4. Corrected sales (v-curve) and final prediction (top), correction coefficient (bottom) (Color figure online)

Fig. 5 .

 5 Fig. 5. 2019 sales predictions for the different methods

Table 1 .

 1 Results on 2019 test data. the week number as an additional input to the other explanatory variables used. Results shown in the following table (

		RMSE MAE MAPE
	SARIMA	67970 56109 31.01
	SARIMAX	73133 62710 34.57
	RANDOM FOREST 63451 55825 31.97
	PROPOSED METH. 56135 43908 23.89
	applied the Random Forest, a popular machine learning algorithm, to our forecasting
	problem by adding

Table


Thomas Milinkovitch
email: t.milinkovitch@hotmail.fr

Stefano Marras

Fabio Antognarelli

Christel Lefrançois

Stéphane Le Floch

Paolo Domenici

Stéphane Le Floch

The effects of hypoxia on aerobic metabolism in oil-contaminated sea bass (Dicentrarchus labrax)

Keywords: oil, hypoxia, Dicentrarchus labrax, metabolic scope, biliary metabolites, active metabolic rate

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

While coastal ecosystems contribute to the economic value and biodiversity of marine environments [START_REF] Costanza R, D'arge | The value of the world's ecosystem services and natural capital[END_REF]Worm et al., 2006), these habitats are under multiple anthropogenic stressors. Hypoxia is one such threatening stressor, as hypoxic areas have exponentially increased, globally, since the 1960s [START_REF] Breitburg | Declining oxygen in the global ocean and coastal waters[END_REF]. This is mainly due to growing population centres that generate excessive nutrient outputs, leading to the eutrophication of local ecosystems [START_REF] Diaz | Spreading dead zones and consequences for marine ecosystems[END_REF]. Nowadays, more than 400 hypoxic ecosystems are recognised and distributed world-wide [START_REF] Selman | Eutrophication and hypoxia in coastal areas : a global assessment of the state of knowledge[END_REF]. In general, major hypoxic events (i.e. dissolved oxygen levels < 0.5 mL/L, [START_REF] Diaz | Marine benthic hypoxia: a review of its ecological effects and the behavioural responses of benthic macrofauna[END_REF] are publicised, as a mass mortality of aquatic organisms often occurs at these locations (e.g. [START_REF] Nakano | Duration rather than frequency of hypoxia causes mass mortality in ark shells (Anadara kagoshimensis)[END_REF][START_REF] Altieri | Tropical dead zones and mass mortalities on coral reefs[END_REF]. However, sub-lethal hypoxic effects, induced by less intense events, are much more frequent [START_REF] Diaz | Marine benthic hypoxia: a review of its ecological effects and the behavioural responses of benthic macrofauna[END_REF]. In this framework, multiple studies have reported on the impact of hypoxia on several fitness proxies, including fish metabolism (e.g. [START_REF] Claireaux | Influence of water temperature and oxygenation on the aerobic metabolic scope of Atlantic cod (Gadus morhua)[END_REF][START_REF] Claireaux | Responses by fishes to environmental hypoxia: integration through Fry's concept of aerobic metabolic scope[END_REF][START_REF] Rosewarne | Measuring maximum and standard metabolic rates using intermittent-flow respirometry: a student laboratory investigation of aerobic metabolic scope and environmental hypoxia in aquatic breathers[END_REF], locomotor performance (e.g. Lefrançois et al., 2005;Lefrançois and Domenici, 2006;Domenici et al., 2007) and fish growth (e.g. [START_REF] Pichavant | Effects of hypoxia on growth and metabolism of juvenile turbot[END_REF][START_REF] Mcnatt | Hypoxia-induced growth rate reduction in two juvenile estuarydependent fishes[END_REF][START_REF] Jeppesen | Effects of Hypoxia on Fish Survival and Oyster Growth in a Highly Eutrophic Estuary[END_REF].

Hypoxia is not the only stressor affecting marine communities. Petrogenic hydrocarbon contamination is prevalent on a worldwide scale, and is derived from different sources such as shipwrecks, oil-spills, operational discharges and offshore exploration [START_REF] Gesamp | Pollution in the open oceans: a review of assessments and related studies[END_REF].

These issues are relevant since new oil extraction techniques ("unconventional techniques") increase the risk of oil released at sea [START_REF] Gesamp | Pollution in the Open Oceans[END_REF][START_REF] Harriss | Arctic Offshore Oil: Great Risks in an Evolving Ocean[END_REF][START_REF] Shearer | Decreasing supplies, increasing risks in oil development in The oxford handbook of energy and society[END_REF]. The effects of petrogenic hydrocarbons have been extensively investigated at several levels of biological organisation, in different taxa [START_REF] Bejarano | Oil spills and their impacts on sand beach invertebrate communities: a literature review[END_REF]. However, marine ecosystems are rarely under pressure from single factors, and biological impacts are more often the result of synergistic environmental and anthropogenic stressors.

The northern Gulf of Mexico is a well-known example of hypoxia and petrogenic hydrocarbon coexistence. Indeed, while organisms are commonly exposed to hypoxia, thanks to nutrient discharge from the Mississippi and Atchafalaya Rivers [START_REF] Diaz | Spreading dead zones and consequences for marine ecosystems[END_REF], they also face petroleum contamination due to oil spills, such as Deep Water Horizon (White et al., 2012) and operational discharges (Wang et al., 2014). In terms of the scientific literature on this subject, evaluating the combined effects of hypoxia and oil is still a new area of research, where more investigation is required. Nonetheless, limited studies are available and have investigated the physiological effects of both stressors in teleost fish [START_REF] Negreiros | Effects of hypoxia and petroleum on the genotoxic and morphological parameters of Hippocampus reidi[END_REF][START_REF] Claireaux | Effects of oil exposure and dispersant use upon environmental adaptation performance and fitness in the european sea bass, Dicentrarchus labrax[END_REF][START_REF] Dasgupta | Hypoxia enhances the toxicity of cor-exit EC9500A and chemically dispersed southern Louisiana sweet crude oil (MC-242) to sheepshead minnow (Cyprinodon variegatus) larvae[END_REF][START_REF] Mager | Combined effects of hypoxia or elevated temperature and Deepwater Horizon crude oil exposure on juvenile mahi-mahi swimming performance[END_REF][START_REF] Milinkovitch | The effect of hypoxia and hydrocarbons on the anti-predator performance of European sea bass (Dicentrarchus labrax)[END_REF]. Among them, some of these studies have challenged organisms with severe levels of hypoxia and/or contamination, sporadically encountered on the field [START_REF] Claireaux | Effects of oil exposure and dispersant use upon environmental adaptation performance and fitness in the european sea bass, Dicentrarchus labrax[END_REF][START_REF] Dasgupta | Hypoxia enhances the toxicity of cor-exit EC9500A and chemically dispersed southern Louisiana sweet crude oil (MC-242) to sheepshead minnow (Cyprinodon variegatus) larvae[END_REF]. Others assessed biological effects at sub-individual levels (e.g. [START_REF] Negreiros | Effects of hypoxia and petroleum on the genotoxic and morphological parameters of Hippocampus reidi[END_REF], consequently, impact predictions on organism fitness are uncertain.

Therefore, we sought to experimentally evaluate the effects of hypoxia on the metabolic scope (MS) of the European Sea bass (Dicentrarchus labrax), previously contaminated with oil. MS was defined as the difference between the active metabolic rate (AMR) and the standard metabolic rate (SMR), and was chosen as a biological parameter of interest, since it estimates the instantaneous energy an organisms allocates to demanding activities (i.e. locomotion, digestion, feeding, etc.), and is therefore deemed an optimal, measurable fitness proxy (Claireaux and Lefrançois, 2007). Moreover, MS is modulated by hypoxic environments [START_REF] Farrell | Defining Hypoxia: An Integrative Synthesis of the Responses of Fish to Hypoxia[END_REF] and oil exposure [START_REF] Mager | Acute Embryonic or Juvenile Exposure to Deepwater Horizon Crude Oil Impairs the Swimming Performance of Mahi-Mahi (Coryphaena hippurus)[END_REF], but no studies have investigated the cumulative effects of both stressors on sea bass metabolism. Sea bass (Dicentrarchus labrax) was chosen as a biological model thanks to several key characteristics; its primary role as a marine predator, its economic and commercial relevance (aquaculture production of 165, 915 tons and capture fisheries of 6,919 tons in 2016; FAO, 2018), and its perceived sensitivity to hypoxia [START_REF] Szlinder-Richert | PAH metabolites in European eels (Anguilla anguilla) as indicators of PAH exposure: Different methodological approaches[END_REF] and hydrocarbons [START_REF] Claireaux | Avoidance threshold to oil water-soluble fraction by a juvenile marine teleost fish[END_REF].

Materials and Methods

Animals

Forty-four juvenile sea bass, aged 1.5-2 years old were used. Fish were provided by Aquanord Ichtus (Gravelines, France). Sea bass were acclimated in 2 m 3 tanks, with an open circulation of sea water under natural photoperiods of six weeks. Dissolved oxygen concentrations, salinity, temperature and pH were monitored daily [93.23 ± 0.77% air saturation, 35.30 ± 0.06‰, 21.28 ± 0.09°C and 7.96 ± 0.02, respectively (mean ± standard error (S.E)]. Fish were fed with Neostart (Le Gouessant, France), and food was withheld 24 h before exposure to contaminants (see below). The average fish length and weight were, 15.06 ± 0.13 cm and 35.27 ± 0.88 g, respectively.

Exposure to oil and total petroleum hydrocarbons (TPHs) measurement

A Brazilian crude oil was used in this study; the PAH concentration was previously described [START_REF] Milinkovitch | The effect of hypoxia and hydrocarbons on the anti-predator performance of European sea bass (Dicentrarchus labrax)[END_REF]. Briefly, the concentration of 21 US-EPA PAHs (parent and alkylated) was 8.264 mg PAHs per g petroleum. The viscosity at 20°C was 35,574 mPa and the pour and flash points were -28°C and -2.5°C, respectively.

The contamination set-up was described by previous studies (Milinkovitch et al., 2011a, b;[START_REF] Claireaux | Effects of oil exposure and dispersant use upon environmental adaptation performance and fitness in the european sea bass, Dicentrarchus labrax[END_REF][START_REF] Milinkovitch | Biomarker modulation associated with marine diesel contamination in the Iceland scallop (Chlamys islandica)[END_REF][START_REF] Dussauze | Innate immunity and antioxidant systems in different tissues of sea bass (Dicentrarchus labrax) exposed to crude oil dispersed mechanically or chemically with Corexit 9500[END_REF][START_REF] Geraudie | First evidence of marine diesel effects on biomarker responses in the Icelandic scallops, Chlamys islandica[END_REF][START_REF] Milinkovitch | The effect of hypoxia and hydrocarbons on the anti-predator performance of European sea bass (Dicentrarchus labrax)[END_REF]. Briefly, a water pump linked to a funnel, continuously sucked up petroleum and/or seawater from the surface and expelled the mixture to the bottom of a 300 L tank. This generated a homogenous dispersion of oil in the water column, after a period of 6 h [START_REF] Dussauze | Innate immunity and antioxidant systems in different tissues of sea bass (Dicentrarchus labrax) exposed to crude oil dispersed mechanically or chemically with Corexit 9500[END_REF].

In this study, 15 g oil (or 15 g seawater for the control treatment) were poured into the funnel and after 6 h, pairs of fish were exposed for 18 h. Every day, two fish were exposed to oil treatments, and two fish were exposed to the control treatment (four fish exposed per day). In total, we carried out eleven days of exposure; 22 fish were exposed to oil treatments, and 22 fish were exposed to control treatments. During the exposure period, temperature, salinity, pH and dissolved oxygen were stable: 21.30°C ± 0.10, 35.06 ± 0.10‰, 8.06 ± 0.01, 99.60 ± 0.40 % air saturation, respectively.

For each exposure, three 125 mL samples were collected from the tank bottom at the beginning (T = 0 h) and at the end of exposure (T = 18 h) to measure TPH concentration.

Chemical assays were conducted according to [START_REF] Fusey | Comparaison de deux méthodes d'évaluation de la biodégradation des hydrocarbures in vitro[END_REF]. Briefly, samples were extracted using 25 mL dichloromethane (Carlo Erba Reactifs, France), thus enabling separation, by decantation of organic and aqueous phases. The organic phase was removed and dried on anhydrous sulphate to avoid residual from the aqueous phase. Then, extracts were analysed at 390 nm, on a UV spectrophotometer (Unicam, France). The limit of quantification was < 1 mg/L.

Respirometry measurements under hypoxic or normoxic conditions

Every day, of the four fish exposed to oil or control treatments, two were individually exposed to hypoxia and two were individually exposed to normoxia. Finally, four groups of 11 fish were generated: 1) fish previously exposed to control treatment, and then normoxic conditions (Ctrl norx group); 2) fish previously exposed to control treatment, and then hypoxic conditions (Ctrl hyp group); 3) fish previously exposed to oil treatment, and then normoxic conditions (Oil norx group) and 4) fish previously exposed to oil treatment, and then hypoxic conditions (Oil hyp group). Practically, in order to expose fish to normoxia or hypoxia, each contaminated fish was immediately transferred with a net (air exposure < 10 s) to a 0.5 L cylindrical static respirometry chamber (Loligo Systems, Denmark), and immersed in a 100 L clean sea water tank (21.50  0.06°C). Before any respirometry measurements, normoxic or hypoxic conditions were established. Normoxic conditions were generated by a constant bubbling of air into the tank. This allowed the oxygen concentration to persist at percentages > 95% air saturation (6.84 mg O 2 /L). In contrast, hypoxic conditions were generated by bubbling nitrogen into the tank, and monitoring levels through an oxygen probe that transferred this information to an oxygen regulator (OXY-REG, Loligo Systems, Denmark), linked to a solenoid valve (Gas control set, Loligo Systems, Denmark). Depending on tank oxygen levels, the regulator allowed the bubbling of nitrogen (i.e. the solenoid valve was open) from a nitrogen gas bottle in a gas equilibration column provided with seawater from the tank (tanks to a water pump EHEIM 600, UK). Deoxygenated water was then sent back into the tank. For these hypoxic conditions, the oxygen regulator was set at 25% air saturation (1.80 mg O 2 /L), resulting in homogenous concentrations between 24-26% air saturation (1.72 mg O 2 /L and 1.87 mg O 2 /L, respectively). Oxygen reduction was performed during 3 h to generate a progressive oxygen decline ( e ran ois and omenici, 2006), to allow fish to acclimate to hypoxic conditions. After 3 h (performed for fish under normoxic and hypoxic conditions), metabolic rates were measured by intermittent flow respirometry: an external pump that initially allowed for the transfer of water from the tank into the respirometer was turned off every 30 min, for a 15 min duration. The respirometer was in a closed state, and oxygen consumption was measured every 15 s with a dipping probe mini sensor (Loligo Systems, Denmark), associated with optodes (Oxy-4 mini; PreSens Precision Sensing GmbH, Regensburg, Germany) and appropriate software (Pre-Sens Oxy 4v2). Oxygen consumption was calculated using linear least-squares regression (excluding the first and last 2 min of each 15 min closed phase) and expressed as mg of O 2 . h -1 .kg -1 . Then, to integrate the log-log allometric relationship between SMR and body mass, data were corrected to a fish body weight of 100 g, using 0.77 as a weight exponent [START_REF] Claireaux | Effect of temperature on maximum swimming speed and cost of transport in juvenile European sea bass (Dicentrarchus labrax)[END_REF]. Fish were left in the respirometer for 24 h, and oxygen consumption was measured every 30 min. SMR was evaluated after 20 h, using the 15 th percentile method [START_REF] Chabot | Environmental hypoxia as a metabolic constraint on fish: the case of Atlantic cod, Gadus morhua[END_REF]. AMR was measured after 24 hwith only 4 h between SMR and AMR measurements, thereby preventing a decrease of PAHs body burden (see 4.1. for more details)-. In order to measure the AMR, fish were removed from the respirometer and placed in a 10 L bucket, where they were chased to exhaustion. Fish were considered exhausted when they did not respond to gentle tail pinching-and then subjected to brief air exposure (30 s) [START_REF] Roche | Finding the best estimates of metabolic rates in a coral reef fish[END_REF]. They were immediately placed in the static respirometer, where oxygen consumption was measured using linear least-squares regression, and data corrected by fish body mass using the log-log allometric relationship between oxygen consumption and body mass.

Finally, MS was calculated for each fish as the difference between AMR and SMR. After AMR measurements, fish were euthanised with eugenol (4-allyl-2-methoxyphenol, Sigma Aldrich, France) at 0.2 mL/L in sea water. The gallbladder was then excised and frozen at -80°C. Bile from the gallbladder was subsequently used to measure PAH metabolites.

Fixed wavelength fluorescence analysis of bile to estimate PAH metabolites

Bile was removed from the gallbladder and diluted in absolute ethanol (Sigma Aldrich) at a ratio of 1:500. Assays were conducted on a spectrofluorimeter (FS5, Edinburgh Instruments, UK), for five fixed wavelength fluorescence (FF) analysis of bile. Naphthalene, pyrene, benzo[a] pyrene-derived metabolites were detected at emission/excitation wavelength spectra of 290 nm/335 nm, 341 nm/383 nm and 380 nm/430 nm, respectively [START_REF] Aas | PAH metabolites in bile, cytochrome P4501A and DNA adducts as environmental risk parameters for chronic oil exposure: a laboratory experiment with Atlantic cod[END_REF].

Fluorene type derived metabolites were detected at 275 nm/328 nm, according to [START_REF] Kopecka-Pilarczyk | Biochemical response in gilthead seabream (Sparus aurata) to in vivo exposure topyrene and fluorene[END_REF]. Phenanthrene type derived metabolites were detected at 250 nm/370 nm according to [START_REF] Szlinder-Richert | PAH metabolites in European eels (Anguilla anguilla) as indicators of PAH exposure: Different methodological approaches[END_REF]. Values were expressed as fluorescence arbitrary units (au).

Statistical analysis

A two-way analysis of variance (ANOVA) was employed to test the effect of "contamination" and "oxygen levels" on SMR, AMR, MS and FF of bile. STATISTICA (version 10, STATSOFT USA) software was used for all analysis, and significance was accepted at P < 0.05. If significant effects (of a single factor, of both factors and/or factor interactions) were highlighted, a Fisher post-hoc test was performed for intergroup comparisons. Student t-tests highlighted statistical differences between TPHs concentrations at T = 0 h and at T = 18 h exposure. Results are expressed as mean ± standard error (S.E).

Results

Hydrocarbon exposure

TPH mean values were 13.6 ± 0.5 mg/L at T = 0 h and 14.8 ± 0.8 mg/L at T = 18 h. We observed no differences in TPH concentrations for the 11 rounds of fish exposure (P > 0.05).

We also observed no differences between mean TPH concentrations at 0 h and 18 h (P > 0.05). Additionally, no TPHs were detected in the water column of the control treatment.

ANOVA results showed that oil significantly affected the relative concentrations of all PAH types (i.e. naphtalene, fluorene, phenanthrene, pyrene, benzo[a]pyrene) (P < 0.01; Figure 1).

Post hoc tests demonstrated significantly (P < 0.01) higher concentrations of metabolites in contaminated fish (Oil norx and Oil hyp), when compared to uncontaminated fish (Ctrl norx and Ctrl hyp). Results showed that hypoxia did not affect the relative concentrations of PAH metabolite types since no significance was highlighted between Oil norx and Oil hyp groups nor between Ctrl norx and Ctrl hyp groups.

Biological effects on metabolic parameters

Our data showed significant effects of oxygen parameters on fish SMR, with mean values for the Ctrl hyp and Oil hyp groups significantly (P  0.01) lower than the Ctrl norx and Oil norx groups (according to post-hoc tests). No effects of contamination, nor interactions of factors on SMR were detected (P  0.05).

The highest AMR was observed for the Ctrl norx group. AMR for the Oil norx group was significantly lower than the Ctrl norx group (P  0.05). Ctrl hyp and Oil hyp groups generated the lowest AMR: both were significantly lower than the Ctrl norx group (P < 0.01 for both comparisons), and both were significantly lower than the Oil norx group (P < 0.01 for both comparisons).

We observed the same patterns for MS; the Ctrl norx group displayed the highest MS; MS was significantly lower for the Oil norx group (P  0.05, when compared to Ctrl norx group); the Ctrl hyp and Oil hyp groups exposed the significantly lowest MS (with P < 0.01 i.e. when each group was compared to Ctrl norx or Oil norx groups).

Discussion

TPHs concentrations and PAHs metabolites

Our results revealed TPH concentration stability in the tank of exposure during the 18 h exposure. In addition, TPH concentrations did not show any statistical differences between replicates, demonstrating repeatability of the contamination protocol. Our results also showed that fish exposed to oil treatment were contaminated with 13.6 ± 0.5 mg/L and 14.8 ± 0.8 mg/L oil, at T = 0 h and T = 18 h, respectively. This order of magnitude agreed with data from [START_REF] Lewis | Oil spill dispersants[END_REF] when modelling oil slick dispersion, and with measurements conducted in the field, during oil spills [START_REF] Lunel | The Braer oil spill: oil fate governed by dispersion[END_REF]. Moreover, it is interesting to note that ships of 400 gross tonnage and above are allowed to discharge oil contaminated water into the sea if concentration does not exceed 15 mg/L of seawater (Convention Marpol, 1973). Thus, the oil concentrations used in this study were appropriate for investigating the impact of oil spills and operational discharge.

In order to investigate the combined effects of oil and hypoxia, we first exposed fish to oil and then to hypoxia for 24h. Although a decrease of hydrocarbon body burden could operate during these 24h, we believe it was not substantial since [START_REF] Dussauze | Innate immunity and antioxidant systems in different tissues of sea bass (Dicentrarchus labrax) exposed to crude oil dispersed mechanically or chemically with Corexit 9500[END_REF] showed that 14 days were required to observe an approximate 50% decrease in total PAHs in the muscle of Dicentrarchus labrax. Be that as it may, we measured relative biliary metabolite concentrations (at the end of the hypoxic/normoxic period) in order to ensure the differential PAH contamination between contaminated and uncontaminated organisms.

Our results showed that all types of metabolites were significantly higher in the bile from fish exposed to oil treatment, when compared to control fish. This validates the uptake by organisms of all types of PAHs studied (i.e. naphthalene, fluorene, phenanthrene, pyrene, benzo[a]pyrene). Bile metabolite incorporation of light PAHs was expected, since these compounds are highly soluble in seawater, and are likely to be incorporated into aquatic organisms by diffusion [START_REF] Neff | Polycyclic aromatic hydrocarbons in the aquatic environment: sources, fate and biological effects[END_REF]. In contrast, the observed incorporation of heavier and the more hydrophobic PAHs, such as benzo[a]pyrene, did not occur by similar processes, since these compounds have a low solubility in seawater [START_REF] Neff | Polycyclic aromatic hydrocarbons in the aquatic environment: sources, fate and biological effects[END_REF]. Most likely, incorporation was due to contact between suspended oil droplets and the gill surface [START_REF] Ramachandran | Influence of salinity and fish species on PAH uptake from dispersed crude oil[END_REF], but we cannot exclude ingestion as a potential mode of incorporation. Benzo[a]pyrene uptake was highly relevant since this PAH has been shown to induce developmental, reproductive and immune toxicity, and is proposed as the most potent carcinogenic PAH molecule [START_REF] Fox | Toxicological Review of Benzo[a]pyrene[END_REF]. Our results showed that fish exposed to hypoxia or normoxia, after oil contamination, had the same levels of PAH metabolites in their bile, suggesting that the hypoxic period after contamination did not reduce or increase metabolization processes, and would therefore not modulate hydrocarbon body burden in fish.

The effects of oil contamination on fish MS under hypoxic conditions

Our physiological data were in line with other studies. For example, [START_REF] Claireaux | Effect of temperature on maximum swimming speed and cost of transport in juvenile European sea bass (Dicentrarchus labrax)[END_REF] measured AMR and SMR values of 408 and 177 mg of O 2 . h -1 .kg -1 , respectively, for juvenile sea bass acclimated to 22°C in control conditions. Although these values were almost half the values observed in our study, they were in the same order of magnitude. This difference could be attributable to different sizes of fish tested and/or to our methodology, since we determined AMR and SMR by static respirometry, whereas [START_REF] Claireaux | Effect of temperature on maximum swimming speed and cost of transport in juvenile European sea bass (Dicentrarchus labrax)[END_REF] used a swimming tunnel.

Our results initially showed an AMR decrease under hypoxic conditions, which was in agreement with previous studies for teleost fish (Claireaux and Lagardère, 1999;Lefrançois and Claireaux, 2003;[START_REF] Dupont-Prinet | An Investigation of Metabolic Prioritization in the European Sea Bass, Dicentrarchus labrax[END_REF]. Indeed, while some compensatory mechanisms exist (e.g. hyperventilation), fish cannot maintain optimal oxygen extraction from seawater, when facing hypoxic conditions [START_REF] Farrell | Encyclopedia of fish physiology: from genome to environment academic press[END_REF]. As a consequence, AMR is reduced and other aerobic performances, such as locomotion, are influenced by such reductions [START_REF] Farrell | Encyclopedia of fish physiology: from genome to environment academic press[END_REF]. Several studies have observed that moderate hypoxia induces AMR reductions in teleost fish: e.g., Lefrançois et al. (2003) and [START_REF] Fitzgibbon | Metabolic scope, swimming performance and the effects of hypoxia in the mulloway, Argyrosomus japonicus (Pisces: Sciaenidae)[END_REF] showed significant AMR reductions when teleost fish (Solea solea and Argyrosomus japonicus) were exposed to 75% air saturation. In our study, we exposed fish to 25% air saturation (1.80 mg O 2 /L at 21.5°C), which induced a substantial significant reduction (67%) in AMR, thereby contributing to a drastic 94% reduction in MS. In accordance with Fry [START_REF] Fry | Effects of the environment on animal activity[END_REF], such a strong reduction in MS prevents fish from conducting discretionary activities (reproduction, locomotion activities, growth, etc.).

Additionally, a significant reduction in SMR (26%) under hypoxia, suggests that fish were clearly below their critical oxygen concentration (CrC), at 1.80 mg O 2 / L. Based on a study by Claireaux and Lagardère (1999) in sea bass, we expected our fish to reach their CrC at 1.80 mg O 2 /L, but not such a decrease below normoxic SMR control values. Based on this reduction of the SMR we can estimate that, this level of hypoxia (25% a.s.), in addition to provoke a drastic reduction in energy for discretionary activities, would also impact on metabolic energy required to maintain mandatory activities, such as respiration and osmoregulation. Thus, these observations suggest that, at this level of hypoxia, survival could be compromised. Based on the literature, it is difficult to estimate if our experimental conditions could impact fish survival on a long or short-term basis. [START_REF] Claireaux | Effects of oil exposure and dispersant use upon environmental adaptation performance and fitness in the european sea bass, Dicentrarchus labrax[END_REF] exposed D. labrax to a drastic decrease in oxygen concentrations (< 20% a.s. after 1 h, and < 5% a.s. after 10 h). Using this approach, all fish were near death (loss of equilibrium) after 12 h. In contrast, Pichavant et al. (2001) observed no death when D. labrax was exposed to 3.2 mg O 2 /L (42.9 % a.s.) over 42 days.

In relation to the impact of physically dispersed oil (normoxic conditions), our results revealed that this contamination induced a 20% decrease in AMR, and a consequent 32% decrease in MS. This data suggested a reduction in the metabolic energy required for discretionary activities. These results are in contrast to a previous study [START_REF] Milinkovitch | Effect of dispersed crude oil exposure upon the aerobic metabolic scope in juvenile golden grey mullet (Liza aurata)[END_REF], that showed no effects of physically dispersed oil on AMR in golden grey mullet (Liza aurata). This discrepancy could be due to different AMR methodologies (static respirometry vs. swimming respirometry) and/or to species specific differences, suggesting that Liza aurata may be more resistant to hydrocarbons than D. labrax. Different oil (experimental contaminant) origins (Brazilian vs. Arabian extraction), and the resulting different physicochemical properties could explain these discrepancies. For instance, different oil viscosities may have impacted oil dispersion characteristics, affecting processes such as oil bioaccumulation or the adherence of oil droplets to gill structures, etc. Although we did not explore the cause of these AMR reductions, we hypothesize that there was a negative impact of oil contamination on the gills, as the first target organ. Indeed, several studies have shown histopathological and functional changes in fish gills exposed to oil [START_REF] Hesni | Gill Histopathological Changes in Milkfish (Chanos chanos) Exposed to Acute Toxicity of Diesel Oil[END_REF][START_REF] Agamy | Sub chronic exposure to crude oil, dispersed oil and dispersant induces histopathological alterations in the gills of the juvenile rabbit fish (Siganus canaliculatus)[END_REF][START_REF] Gagnon | Metabolic Enzyme Activities in Fish Gills as Biomarkers of Exposure to Petroleum Hydrocarbons[END_REF][START_REF] Cohen | Metabolic Responses of Fish Following Exposure to Two Different Oil Spill Remediation Techniques[END_REF]. Since SMR was not impacted by oil, we believe that functional changes may have limited oxygen diffusion when the animal is challenged, i.e. during exercise, but were unimportant enough to alter oxygen diffusion when the animal was resting.

In terms of the impact of contaminated fish under hypoxia (i.e. the interaction of both stressors), our data revealed a significant 67% reduction in AMR, a drastic 95% reduction in MS, and a significant 24% reduction in SMR, when compared to uncontaminated fish under normoxia (i.e. the Ctrl norx group). These observations suggest a reduction in the fish's energy potential to conduct discretionary, as well as mandatory activities. Surprisingly, the impact on these three metabolic parameters was similar, whether fish under hypoxia were contaminated or not, suggesting that at this level of hypoxia, hydrocarbons had no significant effects on fish metabolism, or that the effects of hydrocarbons were masked by this level of hypoxia. This could be due to the fact that when fish respiration is intense (i.e. during AMR in normoxia), hydrocarbon contamination limits this process by limiting oxygen transfer via the gills, whereas when fish respiration is drastically reduced by hypoxia, hydrocarbon contamination may not be a limiting factor. On this basis, it would be interesting to test the effects of hydrocarbons at moderate degrees of hypoxia (i.e. when respiration is not intensively reduced, and hydrocarbons are a limiting factor) to highlight the synergistic or cumulative effects of both stressors in these environments. [START_REF] Mager | Combined effects of hypoxia or elevated temperature and Deepwater Horizon crude oil exposure on juvenile mahi-mahi swimming performance[END_REF] showed a similar effect of hypoxia alone or combined with oil, on the critical swimming speed of Mahimahi, but no effects when fish were exposed to hydrocarbons alone. The degree of hydrocarbon contamination used by [START_REF] Mager | Combined effects of hypoxia or elevated temperature and Deepwater Horizon crude oil exposure on juvenile mahi-mahi swimming performance[END_REF] was lower than our study, which could explain the lack of effect of this contaminant.

In conclusion, our study showed an effect of hypoxia on all metabolic parameters, including MS, suggesting that fitness is impacted at oxygen values ≤ 1.8 mg/L. Additionally, AMR and consequently MS were also decreased due to oil contamination under normoxia. Surprisingly, no synergistic or cumulative effects of both stressors were observed, suggesting that contamination by hydrocarbons impacts fish metabolism when exposed to normoxia, but no additional effects were observed under hypoxia. However, such an interpretation should be cautiously treated, bearing in mind this study investigated the combined effect of hypoxia and hydrocarbons on a single species, at one hypoxia level. Additionally, our experiments dictated that fish were subjected to specific conditions, such as hypoxia duration and the sequential application of stressors (i.e. fish were contaminated and then exposed to hypoxic conditions). Indeed, the combined effects of hypoxia and hydrocarbons might have been different if contamination had occurred during hypoxia.

Therefore, further investigations are necessary to determine the combined effects of both stressors that coexist worldwide. For instance, future studies could investigate how hydrocarbons modify oxygen threshold values below which fish metabolism is impacted.

Furthermore, studies on physiological functions linked to animal fitness (e.g. reproduction, sexual behaviour) would increase our knowledge on the effects of these combined stressors on fish and aquatic organisms. Several other interactions could be investigated, such as temperature with contaminants, especially within the current marine heat wave context.

Acknowledgements

We thank the European Commission, and especially Marie Curie Action for their financial support o the ' ISHY' project, and the individual European postdoctoral research fellowship.

Figure 1. The relative concentrations of polycyclic aromatic hydrocarbon (PAH) metabolites in fish bile, expressed in fluorescence intensity units, for each fish group: 1) Ctrl norx are sea bass previously exposed to a control treatment, and normoxic conditions. 2) Ctrl hyp are fish exposed to a control treatment, and hypoxic conditions. 3) Oil norx are fish exposed to oil, and normoxic conditions. 4) Oil hyp are fish exposed to oil, under hypoxic conditions. Naphthalene, fluorene, phenanthrene, pyrene, benzo[a]pyrene type derived metabolites were measured using fixed wavelength fluorescence (FF) at emission/excitation wavelengths of 290 nm/335 nm, 275 nm/328 nm, 250 nm/370 nm, 341 nm/383 nm and 380 nm/430 nm, respectively. All are represented in histograms A, B, C, D, E. Values represent means ± standard errors of 11 individuals. Different letters above bars indicate significant differences, where P < 0.05. Two-way ANOVA was used for statistical analyses. FI = fluorescence intensity and a.u. = arbitrary units.


J Boissonade

Oscillatory Dynamics Induced in Polyelectrolyte Gels by a Non-Oscillatory Reaction: A Model

Keywords: . 8 2.40.Ck, 82.40.Bj, 82.33.Ln

We develop a general model and the associated numerical algorithm to compute the swelling dynamics of chemo-responsive polyelectrolyte gels immersed in a reactive ionic solution kept at a nonequilibrium stationary state by a permanent feed of fresh reactants. Using an autocatalytic bistable but non-oscillatory reaction, namely, the bromate-sulfite reaction, we predict that a piece of hydrogel that swells/shrinks as a function of pH, can exhibit spontaneous mechanical and chemical oscillations. This constitutes the extension to realistic and experimentally feasible conditions of results previously obtained on a toy model with artificial swelling conditions.

Introduction

The swelling of gels under the action of external physical or chemical parameters is raising a growing interest in view of multiple potential applications in sensors, actuators, biomimetic systems or drug delivery. Among swelling agents able to induce swelling, the chemical composition of the solvent has received particular attention [START_REF]Nonlinear Dynamics related to Polymeric Systems[END_REF][2][START_REF]Nonlinear Dynamics in Polymeric Systems[END_REF][4][START_REF]Chemomechanical Instabilities in Responsive Material[END_REF]. In this article, we consider the coupling of chemical reactions with chemo-responsive gels as a source of spontaneous mechanical oscillations. Different strategies to get such mechanical oscillations have been used. In all cases, a small piece of hydrogel that swells/shrinks as a function of the composition of the solvent is plunged in a solution of chemicals in large excess or permanently fed with fresh reactants. These reactants diffuse into the gel where the reaction is going on and drive the swelling process. More extensive reviews than the following short introduction will be found in [START_REF]Chemomechanical Instabilities in Responsive Material[END_REF]. First, one can immerse a piece of gel in an oscillating reaction carried out in a continuous stirred tank reactor (CSTR) or a large reservoir. The periodic changes of concentrations induce periodic swelling/shinking cycles of the gel [START_REF] Yoshida | [END_REF]7]. The process was simulated with toy models by Villain et al. [8,9]. Alternatively, one can attach to the gel a catalyst necessary for the reaction to oscillate by grafting it on the polymer network. This has been done for the Belousov-Zhabotinskii reaction (B.Z.) by Yoshida and coworkers [10][11][12] and a model has been proposed by Yashin and Balasz [13]. In this case, only the gel and its contents are oscillating since the required conditions are only fulfilled within it. However, in both cases, the source of the oscillations is a purely chemical instability. Mechanical oscillations are slaved to chemical oscillations that would also occur in an inert gel (apart from slight changes due to the dilution of the catalyst during swelling). Unfor-tunately, oscillating reactions remain rare and the domain of parameters where oscillations occur is often very small.

Other strategies based on bistability and on the associated hysteresis, a common phenomenon in nonlinear physics, have been proposed. The basic principle is to create a feedback that makes both states slightly unstable so that the system permanently switches between them, following periodically the hysteresis cycle. Siegel et al. use the hysteresis of the volumic transition in a membrane made of a gel responsive to pH to modulate its porosity and the coupling through this membrane of two compartments where an autocatalytic reaction takes place, leading to chemical oscillations [14][15][16]. Instead, we proposed to take advantage of chemical bistability and to use the volume changes caused by swelling/shrinking as a feedback [17][18][19]. This was partly supported by a series of experiments performed with the chlorite-tetrathionate (CT) reaction in a poly(NiPAA-co-AA) gel [20][21][22]. We shall now focus on this approach.

Although our former work [17][18][19][20][21][22] have clearly demonstrated the validity of chemical spatial bistability as a potential source of oscillations, there were still two sorts of problems to be solved to get a coherent whole. First, even in a non-responsive gel, the CT reaction exhibits, in addition to spatial bistability, complex dynamics [23] which perturbs the interpretations. Moreover, although a description in terms of effective diffusion coefficients is sufficient to understand the essential dynamics, the kinetic model breaks down when the motion of ionic charges is accounted for. The choice of a more apropriate reaction is advisable. The second problem is also fundamental. In our former toy models, the gel swelling properties were related to the chemical composition of the solvent in an heuristic way through an arbitrary dependence of the Flory parameter χ with one of the concentrations. More satisfactory models should rely on a more pertinent description of the swelling source. In pH sensitive polyelectrolytes, the most commonly used responsive hydrogels, the swelling is induced by the osmotic pressure due to the partial ionisation of the polymer network. The purpose of this article is to demonstrate that occurence of oscillations that originate in pure chemomechanical instabilities based on spatial chemical bistability can be predicted on the basis of reasonable and more realistic models. We shall use the Bromate-Sulfite reaction (BS) for which a good kinetic model is available [25,26] and we shall develop a model of swelling for a polyelectrolyte in presence of the reaction. Beforehand, we shall briefly summarize the principles of the chemo-mechanical instability. For more extensive developments on spatial bistability see [23,24] and on the chemo-mechanical instability itself, see [18,19,23].

Many autocatalytic reactions carried on in a CSTR exhibit bistability as a function of the residence time t res [START_REF] Epstein | An Introduction to Nonlinear Chemical Dynamics[END_REF] . At short t res , the extent of reaction ξ is low. The reaction time is longer than the time to refill the reactor with fresh reactants and the composition remains close to the concentrations in the input flow. On the opposite, for a long t res , ξ is large and the reaction is almost completed. At intermediate values of t res , there are two stable states, one with low ξ (the "flow" state F) and the other with large ξ (the "thermodynamic" state T) selected according with the previous history of the system. If a piece of a chemically inert gel, plunged in the same almost unreacted medium -most often a CSTR kept in the F state -a similar phenomenon called "spatial bistability" occurs as a function of the typical size l of the gel. This parameter, which determines the time to transport by diffusion the reactants to the deep core of the gel, plays a role analog to t res . If l is small, exchanges with the environment are faster than reaction rates. The whole gel remains in a state close to the reservoir (low ξ) and is again referred as a flow state (F). If l is large, the reaction becomes dominant at some distance of the boundary, so that the deep core of the gel is at large ξ, except for a boundary layer that insures continuity of concentrations with the reservoir. This defines the "mixed state" M. At intermediate sizes l inf < l < l sup , both states can be stable, which defines the spatial bistability. A number of reactions have already been shown to exhibit this phenomenon [23,26,[START_REF] Blanchedeau | [END_REF][29][30]. Obviously, other parameters such as a concentration in the feeding medium can be used to define the bistability domain, but l is the most relevant to our goal. If one could slowly and continuously increase l, starting from l < l inf , one could follow the F state branch until the system switches to the M state at l = l sup . If, now, we reverse this change, l decreases and the system follows the M branch down to l = l inf , where it switches back to the F state. In this way, one follows the whole hysteresis cycle associated to the spatial bistability. If we replace the inert gel by a chemo-responsive gel, the changes of size l are spontaneously driven by the chemical state. For instance, in the BS reaction that will be studied further, the F state is alkaline and tends to make the polyacid gel to swell, whereas the M state is mainly acid and tends to make the gel to shrink. If l is appropriate and the swelling/shrinking process is slow enough in regard to the reaction-difusion process, one can start on state F, swell up to l = l sup , switch to state M, shrink down to l = l inf , and switch back to state F. The process is repeated indefinitely, leading to periodic oscillations both of the volume of the gel and of the concentrations within the gel. One does not need an oscillating reaction, only a much more common bistable one, to produce such a chemomechanical instability. Some models have been proposed to predict the dynamics of swelling in neutral gels without reaction [31][32][33][34][35][36][37] or with reaction [9,13,18,19]. In the second case, the chemical sensitivity of the gel is descibed by an assumed and somewhat arbitrary dependence of the Flory parameter χ, which gathers all the energetic interactions of monomer and solvent molecules, on the concentration of a chemical species. However, most gels that exhibit large volume variations are polyelectrolytes that mainly swell under the action of the ionic osmotic pressure [START_REF] Flory | Principles of Polymer Chemisty[END_REF][START_REF] Ricka | [END_REF][START_REF] Shröder | Physical Properties of Polymeric Gels[END_REF]. Only a few dynamical complex models have been developed to simulate the swelling of polyelectrolytes in presence of simple ionic reactions [START_REF] Achilleos | [END_REF]42]. With complicated reactions implying both ions and neutral species, the situation is still more complex in regard to the number of involved primary and secondary phenomena with unknown parameters. Thus, to work out a model that could be applied to a realistic reaction, namely, the BS reaction, we had to resort to some simplifications.

In this article, we consider a piece of gel immersed in a CSTR kept in the F state by a permanent feed of fresh reactants. The equations are general but will be eventually applied to the BS reaction. We limit ourselves to a one dimensional (1-D) system. A gel swelling in a capillary is not very realistic from an experimental point of view but a good approximation to such a 1-D system could be obtained by grafting a flat piece of gel on the bottom of the reactor, provided that the thickness l is much smaller than the size in the transverse directions for the effects of the deformations in these directions to be limited to the edges and to be neglected. In Section 2, we introduce a simple model for the mechanism of swelling of the polyelectrolyte in presence of a reaction. In Section 3, we establish generalized reaction-diffusion equations within such a gel. In Section 4, are discussed some aspects of the spatial bistability of the BS reaction. In section 5, we gather all the elements of the model and perform the numerical simulation which demonstrates that chemomechanical oscillations of the type described above can actually be predicted. In an inert gel, the model previously developed for the kinetics and transport in the BS reaction is in good agreement with experimental data [26]. Its extension to the swelling gel remains more likely valid, but modeling the swelling process itself needs a simplified description in which only the main contributions are accounted for. Thus, it must be clear that we do not claim we give an almost exact quantitative modeling of the dynamics, but we believe that this approach captures the main phenomena, the orders of magnitude, and the general trends. It should be considered as a first guide for experimentalists.

The swelling model

There are three type of constituents in the gel: the solvent (water), the polymer, and the solutes (the chemicals).

The concentrations c i of the chemicals are low so their volume fraction is assumed to be negligible in regard to the volume fraction φ of the polymer and the volume fraction 1 -φ of the solvent. The polymer is assumed to be a polyacid. The concentration of ionisable sites in the polymer in the absence of solvent would be c a0 , so that their concentration in the gel is given by c a = φ c a0 . In the referential of the laboratory, the dynamics at a given point in the gel are characterized by the velocity of the polymer matrix v P . The permeation flow of the solvent through the polymer matrix is associated to the relative velocity δv = v P -v S where v S is the velocity of the solvent in the same referential. In common hydrogels, the total net-work+solvent system satisfies the incompressibility condition. Then, in systems with stagnation points or fixed walls (like here), the equation for mass balance expressed in terms of volume fractions is given by the conservation law:

φv P + (1 -φ)v S = 0 (1)
which allows to express v S and δv as a function of v P .

The local force density acting on the polymer matrix is given by f = ∇ • σ where σ is the stress tensor. These forces result from the form of the free energy associated to the mixing of the polymer molecules with the solvent molecules, from the elastic forces exerted by the network when deformations are applied, and from those associated to the partial ionisation of the network. The stress tensor can be splitted in two parts that correspond respectively to an isotropic contribution and an anisotropic contribution in the form:

σ ij = -δ ij Π isotropic + σ (noniso) ij nonisotropic (2)
where Π can be considered as a pressure and is commonly refered to as the "osmotic pressure". It can be decomposed in three terms corresponding to the three sources of forces:

Π = Π mix + Π elas + Π ion (3)
For a unidimensional deformation, there is only one coordinate and the anisotropic part of the stress tensor, which is only due to the elastic forces, can be included in the elastic contribution Π elas to the osmotic pressure. Then, the system is at equilibrium when Π = 0. Out of equilibrium, the osmotic forces are counterbalanced by the friction forces exerted by the solvent on the network and proportional to their relative velocity δv:

Friction forces = Osmotic forces ζ(φ)δv = -∇Π (4)
where ζ(φ) is a friction coefficient which increases with φ. The exact shape of this function is not well known. Many authors use a dependence in φ 3/2 on the basis of theoretical arguments only valid for semidilute polymers.

At the low polymer densities on which we focus, we have preferred the Ogston model based on flow analysis in fiber networks [START_REF] Ogston | Proc. R. Soc. Lond. A[END_REF], as proposed in [34]:

ζ(φ) = RT V S 1 D 0 φ e η √ φ (5)
where η is of the order of the ratio of the radius of the polymer chains to the size of the solvent molecule and V S the molar volume of the solvent. In the original theory, D 0 should be the autodiffusion coefficient of this solvent but it was found in swelling experiments that the actual value should be almost two orders of magnitude larger [34]. Thus, we shall consider D 0 as an expandable parameter in the range 0.01 mm 2 /s ≤ D 0 ≤ 0.1 mm 2 /s. The mixing term is written in agreement with the classical Flory-Huggins theory [START_REF] Flory | Principles of Polymer Chemisty[END_REF][START_REF] Doi | Introduction to Polymer Physics[END_REF]:

Π mix = - RT V S [φ + log(1 -φ) + χφ 2 ] ( 6
)
The two first terms depend on the entropy of mixing, the third one depends on the mutual energetic interactions between the different molecules (solvent and monomers).

A more rigorous theory would need to introduce virial coefficients specific to each gel [START_REF] Grossberg | Statistical Physics of Macromolecules[END_REF]. Large swelling generally occurs when χ is slightly larger than 0.5. There has been a great number of debates about the form of the elastic terms without definite and universal answer. At our level of description and in the absence of precise experimental data -to be determined for each composition of gel -we shall follow the simple phenomenologic approach of Barrière and Leibler [35] based on the work of Bastide and Candau [START_REF] Bastide | Structure of gels investigated by static scattering techniques[END_REF]. The authors assert that both the free energy by unit volume and the elastic modulus G(φ) scale like (φ/ψ) n . To fit at best with other theories and recover Flory's expressions [START_REF] Flory | Principles of Polymer Chemisty[END_REF], we use n = 1/3 and ψ = φ 0 where the reference state of volume fraction φ 0 is supposed to be in a relaxed state (normally a condensed state before swelling). During swelling, the network experiences deformations and a point i of initial coordinate X i is moved to coordinate x i at time t. In the following, we consider the relaxed state as an initial undeformed state. To remain unidimensional during swelling, the gel has to be constrained by external walls or grafting. This implies that, contrary to a gel that swells freely in its solvent, the final equilibrium is not an isotropic state. The local deformation at point i and time t is measured in a 1-D system by dx i /dX i . The unidirectional transformation dX i → dx i can be decomposed in a virtual isotropic swelling dX i → dx ′ i followed by an anisotropic stretching dx ′ i → dx i imposed by the constraints at constant volume with deformation λ = dx ′ i /dx i . Thus

φ 0 φ = dx i dX i = dx i dx ′ i dx ′ i dX i = λ φ 0 φ ) 1 3 ⇒ λ = φ 0 φ ) 2 3 (7)
On the basis of the above assertions, the following expression can be derived [35]:

Π elas = -K net RT V S φ φ 0 1 3 1 + C λ λ 2 - 1 λ (8)
where K net is a constant which depends on the network properties and C λ is a constant of order unity proportional to the shear modulus G(φ). The terms that contain λ = (φ 0 /φ) 2 3 represent the nonisotropic contribution to the osmotic pressure.

To give an expression for Π ion , we shall first consider the equilibrium. In accordance with many authors [START_REF] Shröder | Physical Properties of Polymeric Gels[END_REF], we neglect, in this simple approach, all the electrostatic interactions between the ions in the network to keep only the entropy terms that should be dominant. We closely follow the approach of Ricka and Tanaka [START_REF] Ricka | [END_REF]. In the gel, the ionisable monomers HA, the concentration of which is c a0 φ in the absence of dissociation, are actually partly ionised into H + and A -according to:

HA ⇆ H + + A - (9)
with the equilibrium constant:

K g = [H + ][A -] [HA] (10)
From equation (10), one gets:

c a = K g c a0 φ K g + [H + ] (11)
Whereas other ions are mobile, the ions A -are attached to the network. According to the classical theory of the Donnan equilibrium, this creates an excess of mobile ions in the gel in regard to their concentrations in the surrounding medium (here, the CSTR contents). This creates in turn a pressure difference, which, for small concentrations is analog to the pressure created by a perfect gaz:

Π ion = RT i (c i -c ′ i ) ( 12
)
where the sum was taken on the mobile ions, the c ′ i 's are the concentration of these ions in the CSTR and the affinities have been replaced by the concentrations. To cross the boundary, the ions have to go through a potential difference δU in agreement with the Boltzman law:

c i c ′ i = exp - z i e δU k B T = K zi ( 13
)
where e is the charge of the electron, z i is the charge number of ion i and K is the Donnan ratio, which is common to all mobile ions. Given the c ′ i 's, and using equations [START_REF] Yoshida | [END_REF]8,11,12,13), one can compute K and φ by solving numerically the system of two nonlinear equations which respectively express electroneutrality in the gel and mechanical equilibrium:

i c i z i + z a c a = 0 (electroneutrality) (14)
Π = 0 (mechanical equilibrium) (15) where the charge number of A -is z a =-1. The concentrations of non ionic species are continuous at the boundary.

To solve the dynamical equations, we assume that the system it at local Donnan equilibrium at the CSTR boundary so that, at this point, Π ion , φ and K can be computed as explained above. Consequently, all the c a and c i 's concentrations are obtained from equations ( 11) and ( 13). This defines completely the boundary condition. Inside the gel, the generalisation of the calculation of Π ion is not straighforward, because the reaction creates or destroys not only ions but also neutral molecules. One can assume that the ionic osmotic pressure is proportional to the degree of ionisation of the polymer which can be written in the form:

Π ion = K ion c a (16)
To remain consistent, at the CSTR boundary, Π ion computed from equation ( 16) must be equal to the value already computed from equation (12). Thus, the value of K ion is obtained from equation ( 16) where c a and Π ion at the boundary are respectively obtained from equations ( 11) and (12).

If the distribution of concentrations and φ at time t is known, it is now possible to compute Π at each point inside the gel. Then, from equations ( 1) and ( 4), one gets the velocity of the polymer

v P = - (1 -φ) ζ(φ) ∇Π (17)
and the velocity of the solvent

v S = - φ (1 -φ) v P (18)
Taking into account the equation of conservation for the polymer ∂φ ∂t

+ ∇ • (v P φ) = 0 ( 19
)
one gets the dynamical equation for φ

∂φ ∂t = ∇ φ(1 -φ) ζ(φ) ∇Π (20)
As equation (17), equation (20) describes the motion of the polymer matrix through the evolution of the gel density. To solve our problem, one needs now to couple equation (20) to the equations that describe the reaction and the transport of the chemical species in the gel.

The transport model

Inside the gel, the species are transported by three processes, namely, migration due to the internal electric field, diffusion, and convective motions due to swelling. Our main hypothesis is that the solutes, indexed by i or k, are convected like the solvent (i.e. at velocity v S ), whereas the anion A -, indexed by a, is attached to the network and is obviously convected like the polymer (i.e. at velocity v P ).

The transport model is a straightforward generalisation of the transport equations for electrolytes as presented by Newman [START_REF] Newman | Electrochemical Systems[END_REF]. The fluxes are given by the equations

N i = Nmi -F z i u i ∇U -D i ∇c i + c i v S (21) N a = c a v P ( 22
)
where F is the Faraday constant, u i is the mobility of species i, D i its diffusion coefficient, and U is now the diffusion potentiel created by the distribution of charges when the diffusion coefficients of the ions are not all equal [START_REF] Newman | Electrochemical Systems[END_REF]. In the expression of N i , the first term N mi corresponds to migration, the second term to diffusion, and the late term to convection. The current density is

j = F ( i z i N i + z a N a ) (23)
Taking into account equations ( 21) and ( 22) and the electroneutrality condition

i z i c i + z a c a = 0 (24)
one gets

j = -κ∇U -F i z i D i ∇c i + F δvz a c a (25)
where the conductivity κ is defined by κ = F 2 i z 2 i u i c i . There is no permanent current, so that, writing j = 0 in equation (25), one can eliminate the potential U from equation (21). We define the transference numbers

t i = z 2 i D i c i i z 2 k D k c k (26)
where the diffusion coefficients replace the mobilities since, according to the Nernst-Einstein relation one has D i = RT u i . The migration flux N mi becomes:

N mi = t i z i k z k D k ∇c k -z a c a δv (27)
One can also eliminate one ion, indexed by n, from this expression by rewriting the electroneutrality condition under the form:

-z n c n = k =n z k c k + z a c a (28)
It is advisable to use a non reactive ion with a well known diffusion coefficient D n . This is generally Na + or K + which are often the cation of input reactive salts but which, most of the time, remain spectators of the reaction, being only involved in the transport by their electric charge. The final expression for N mi is:

N mi = t i z i   k =n z k (D k -D n )∇c k -D n z a ∇c a -z a c a δv   (29)
The diffusion coefficients of the solutes depend on the volume fraction φ of the polymer. The diffusion coefficient is proportional to the ratio ǫ/τ of the permitivity to the tortuosity. Since we shall only consider small polymer densities φ ≪ 1, we assume that both ǫ and τ depend linearly on φ. The available volume fraction for the solutes is 1-φ, so that we can write ǫ = 1 -φ. For the tortuosity, one has τ = 1 + αφ where α is a coefficient of order unity (3/2 for a periodic network of spheres [START_REF] Cussler | Diffusion: Mass Transfer in Fluid Systems[END_REF]). For diffusion, the network can be seen as a random network of fibers. Extrapolating the numerical simulations of Tomadakis and Sotirchos [START_REF] Tomadakis | [END_REF] to φ → 0, one find α ≈ 1. Thus, if D 0i is the diffusion coefficient of species i in pure solvent, we get:

D i (φ) ≈ D 0i 1 -φ 1 + φ ≈ D 0i (1 -2φ) (30)
In the computations, we shall use

D i (φ) = D 0i (1 -2φ).
The evolution of the concentrations are given by the conservation laws

∂c i ∂t = -∇ • N i + R i ( 31
)
∂c a ∂t = -∇ • N a + R a ( 32
)
where R i and R a are the reaction terms ruled by the chemical kinetics and N i and N a are obtained from equations ( 21), (22), and ( 29). These transport equations are general. In the next section, we introduce the kinetic terms R i and R a for the BS reaction.

The reaction model 4.1 Kinetic equations for the BS reaction

The Bromate-Sulfite model is a reaction which is autocatalytic with H + and has been studied previously both in CSTR [25] and in the context of spatial bistability [26]. The reader is invited to refer to [26] for details on the kinetic model and the experimental procedures. We shall here report the essentials to understand the next section. The gel is in contact with the contents of a CSTR fed with solutions of bromate, sulfite and sulfuric acid at a residence time τ =500 s. The volume of the CSTR is assumed to be large enough in order that the reaction inside is not significantly influenced by the gel contents. The kinetics equations, kinetic constants, and diffusion coefficients are the same as those used in [26] except for an additional reaction to account for the dissociation of the gel according to equation (9). The pK of the polyacid has been fixed to 5.5 which should corresponds to the pK of a gel of poly(N-isopropylacrylamide-co-acrylic acid) (NiPAAMco-AA) which is commonly used in swelling and chemomechanics experiments. The kinetics are described by a series of balance equations completed by fast equilibria

BrO - 3 +3HSO - 3 → 3SO 2- 4 +Br -+3H + (R1) BrO - 3 +3H 2 SO 3 → 3SO 2- 4 +Br -+6H + (R2) BrO - 3 +6H 2 SO 3 → 3S 2 O 2- 6 +Br -+6H + +3H 2 O (R3) SO 2- 3 +H + ⇆ HSO - 3 (R4) HSO - 3 +H + ⇆ H 2 SO 3 (R5) SO 2- 4 +H + ⇆ HSO - 4 (R6) HSO - 4 +H + ⇆ H 2 SO 4 (R7) H + + OH -⇆ H 2 O (R8) H + + A -⇆ HA (R9)
with the corresponding rate laws

r 1 = k 1 [HSO - 3 ][BrO - 3 ] r 2 = k 2 [H 2 SO 3 ][BrO - 3 ] r 3 = k 3 [H 2 SO 3 ][BrO - 3 ] r 4 = k 4 [SO 2- 3 ][H + ] -k -4 [HSO - 3 ] r 5 = k 5 [HSO - 3 ][H + ] -k -5 [H 2 SO 3 ] ( 33
)
r 6 = k 6 [SO 2- 4 ][H + ] -k -6 [HSO - 4 ] r 7 = k 7 [HSO - 4 ][H + ] -k -7 [H 2 SO 4 ] r 8 = k 8 [H + ][OH -] -K w r 9 = k 9 [H + ][A -] -k -9 [HA]
The constants and the diffusion coefficients in pure water used in the numerical simulations are respectively given in Table 1 and Table 2. When they were unavailable in the literature, the diffusion coefficients were fixed to reasonable values and are marked with an asterisk. Results on spatial bistability in a non responsive and chemically inert gel have been reported in [26]. The introduction of Table 1. Kinetic constants used in the numerical simulations reaction (R9) could modify the results. In order to well distinguish effects of the sole chemistry from those of chemomechanical instabilities, and also check that the parameters requirements for the emergence such intabilities are still fulfilled, we have repeated some of the previous simulations when a non diffusive complexing agent such as A -is present inside the non-responsive gel, but not in the CSTR. In these preliminary computations, the purpose is to check the sole possible influence of the presence of the complex HA. Thus, as in [26], the volume fraction of the gel is neglected. The concentration [HA] 0 before dissociation was fixed to 5×10 -2 M, that will be the typical mean value for c a0 φ that will be obtained for the mechanically responsive polyelectrolytes in the next section.

k1 = 0.0653 M -1 s -1 k2 = 18 M -1 s -1 k3 = 0.7 M -1 s -1 k4 = 5 × 10 10 M -1 s -1 ; k-4 = 3 × 10 3 s -1 k5 = 2 × 10 8 M -1 s -1 ; k-5 = 3.4 × 10 6 s -1 k6 = 1. × 10 10 M -1 s -1 ; k-6 = 1. × 10 3 s -1 k7 = 1 × 10 11 s -1 ; k-7 = 1.148 × 10 9 M -1 s -1 k8 = 1.4 × 10 11 M -1 s -1 ; Kw = 10 -14 × k8 k9 = 1 × 10 10 M -1 s -1 ; k-9 = 10 -5.5 × k8 M -1 s -1

Spatial bistability of the BS reaction in presence of a complexing agent

The concentrations in the input flow, indexed by 0, are always chosen in order that the CSTR remains in the flow state, i.e. slighly basic. The F and M states in the gel are characterized by the pH at the bottom of the gel: basic or quasi neutral for state F, acid (pH∼ 3) for state M. The concentrations at the CSTR/gel boundary (x = l) are derived in agreement with the rules of Donnan equilibrium. The concentrations in the CSTR being first computed, the Donnan ratio K is computed according to the procedure described in Section 2 to determine the fixed concentration within the gel at the CSTR boundary. Since φ = 0, only equation 14 has to be solved with

c a = K g [HA] 0 /(K g + [H + ]
) and the c i 's are given by equation 13. No flux boundary conditions are applied at the opposite wall (x = 0). We found that the presence of the complexing agent does not preclude spatial bistability. In Figure 1, is shown a typical bistability diagram in the plane ([H 2 SO 4 ] 0 , l) where l is the size of the 1-D system. For a given set of input concentrations, the range of sizes

[BrO - 3 ]0=25 mM, [SO 2- 3 ]0=60 mM, [H2SO 2- 4 ]0=6.5 mM, l=1mm)
l for which the system is bistable is an essential parameter for the emergence of chemomechanical oscillations. We found that these limits of the bistability region are almost independent of [HA] 0 in the range of concentrations we shall consider in the following section.

For information, we also report in Figure 2 the concentration profiles of H + and c a in states F and M for a set of data within the bistability range. The system was previously found to exhibit marginal oscillations in a very narrow domain of parameters for larges values of [BrO - 3 ] 0 . However, it is known that, for most oscillatory reactions, complexation of the catalyst generally kills oscillations whereas bistability is preserved. This is actually the case. We have checked that these marginal oscillations are totally killed in the polyelectrolyte gel and that the oscillations that will be evidenced in the next section are of pure chemo-mechanical origin.

Chemo-mechanical oscillations

We now couple the chemical reactions with gel swelling, solving simultaneously equations ( 20), (31), (32) where the R i 's and R a are given by equations (34) and the expressions of ∇Π, N i , and N a are computed as explained in the previous sections.

Numerical techniques

To solve these equations, it is appropriate to switch to a Lagrangian representation on a fixed grid with a uniform volume fraction φ 0 , size l 0 , and a constant spatial stepsize δX i . In this 1-D system, the transformation from eulerian coordinates x to Lagrangian coordinates X is operated by the simple transformation

x -→ X =⇒ ∇ x -→ φ φ 0 ∇ X (34)
in the equations. At each time, one can recover the value of x i corresponding to the point of coordinate X i in the fixed grid from the equation

x i = Xi 0 φ 0 φ(X) dX (35)
to get the values of c i (x i ) and φ(x i ) in eulerian coordinates from c i (X i ) and φ(X i ). The integration on the Lagragian fixed grid is performed with finite differences and a method of lines based on a stiff fourth order Rosenbrok temporal integrator [50].

Parameters choice

We limit ourselves in this presentation to the complete determination of the oscillatory domain in a ([H 2 SO 4 ] 0 ,l 0 ) plane of the phase diagram for a unique representative set of the other parameters. Chemical parameters are archetypal of those that were used in spatial bistability studies of the BS reaction. As for computational results reported in Figure 1, the bromate and sulfite concentrations in the input flow are respectively fixed to [BrO - 3 ] 0 =25 mM and [SO 2- 3 ] 0 =60 mM. Only the size of the gel (through of l 0 ) and [H 2 SO 2- 4 ] 0 are systematically varied.

The choice of gel parameters is more problematic and is only guided by the necessity to ensure that both the values of these parameters and the swelling amount are similar to those encountered in typical experiments. The reference volume fraction φ 0 is assumed to be the value of φ in a fully contracted gel (acid medium). Since the gel must be diluted (typically a few percent in the experiments), we take φ 0 = 0.05. To account for the shearing effects, we use equation 8 with C λ =1 as in [35]. To produce a significant swelling amount without a full volume transition, the Flory parameter is fixed to χ=0.515, i.e. a value slightly larger than 1/2, and the corresponding value of K net is obtained in the relaxed state (no stress) from

K net = -[φ 0 + log(1 -φ 0 ) + χφ 2 0 ]φ -1 3 0 . ( 36
)
As in section 4, the pK a of the gel is fixed to 5.5, the value generally retained for the poly-Nisopropylamide-co-acrylic acid (NiPAAM-co-AA) gels. In the simulations reported in the next paragraph, the coefficients of equation ( 5) are fixed to D 0 =0.03mm 2 /s and η = 5. In these computations, the size of the gel is changing during time. As a size reference for a given system, we use l 0 , the size in the fully contracted state, i.e. at volume fraction φ 0 . Since there is always a part of the gel which is in a swollen state, at least in the region close to CSTR/gel boundary, the real size l(t) at a given instant t is always larger than this value.

Simulation results

The domain where oscillations are actually observed is mapped out as a black region in the ([H 2 SO 4 ] 0 , l 0 ) parameter plane (Figure 3). In the definition of this region, the ordinate l must be understood as the reference length l 0 . The actual range of sizes covered during the oscillations for three representative sets of parameters (white points) is represented as vertical bars. In reference to the heuristic approach, the limits of spatial bistability already computed in section 4 for the same parameters are also replicated from Figure 1. One must be aware that they are given for information but that they should be considered with caution. They were obtained for a stationary system and a fixed value of [HA] 0 . In the present situation, the polymer volume fraction is not uniform. It has some influence on the diffusive transport, on the distribution of charges and additional convective motions are present. Moreover, the swelling time is not infinite in regard to the reaction/diffusion characteristic time. However, one can see that the minima and maxima of size during an oscillation are close the bistability limits, in agreement with the heuristic theory. As expected, the oscillations are periodic. The period T and the amplitude ∆l of oscillations, for the three representative points in Figure 3 are reported in Table 3. To get a better idea of the effective size of the system, we have added l m , the median value between the maxima and minima and the relative amplitude ∆l/l m (expressed in percent). For illustration, both the oscillations of size and those of the concentration [H + ] at the fixed wall are shown in Figure 4 for point (2). Coupling the gel and the reaction not only induces mechanical pulsations but also gives rise to large pH oscillations in the deep core of the gel, whereas the CSTR/gel interface keeps fixed concentrations.

When the size increases, the period changes from a few hours to several days, so that using sizes l m > 1 mm would not be realistic from an experimental point of view. Close to the limits of the bistability region in a non-responsive gel, i.e. at the largest [H 2 SO 4 ] 0 values (see point 1) for which oscillations are observed, the amplitude of oscillations decreases as does the concentration gap between the two states. In this region, whereas the swelling time decreases with the system size, reaction dynamics exhibit some slowing down at the approach of the critical point so that the two characteristic times are getting closer and the oscillations present less relaxational features as can be seen on Figure 5 where we compare the shape of an oscillation at point (1) to the shape of an oscillation at point (3).

Discussion and conclusions

Although our analysis is applied to a 1-D system, one expects that qualitative behavior and orders of magnitude should be similar in other geometries since the increase in size ratio that could result from constraining the system to swell in a particular direction is somewhat counterbalanced by a larger shear stress that reinforces the elastic forces which limit swelling. Actually, preliminary experiments performed in our laboratory in conical gels evidence mechanical oscillations with a period of a few hours [51]. More elaborate experimental systems are presently developed to get confirmation and reliable data for a more systematic analysis. Examination of Figure 3 provides important informations for experimentalists. An important result is that, for [H 2 SO 4 ] 0 lower than a value close to the critical bistability limit, there is always a small range of sizes in which the gel exhibits mechanical oscillations. Another point is that, even close to the critical bistability limit where small sizes allow for the shortest periods, the most suitable for experiments, the relative amplitude of oscillations remains large enough for observations. Although the model seems able to provide good qualitative predictions and general trends, there are some important limitations and approximation levels that should be considered in the future. To conclude, let us briefly summarize the most important assumptions that could limit the validity of the model or its extensions to other problems. We have neglected the role of electrostatic interactions within the polymer chains. Their role is controversial, although most authors consider the entropic effect of the excess of free charges to be dominant. However, for strong densities of charges, counterions condensation [52] could significantly modify swelling predictions. We have, as usual, also replaced activities by concentrations but the domain of validity of this approximation is smaller for ionic solutions than for neutral species. We also assumed than the solutes are convected as the solvent and that their diffusion is only modified by steric hindrance by means of equation (30), neglecting all cross terms that would express other molecular interactions. By using a local Donnan equilibrium at the CSTR/gel boundary, we also implicitly assumed that the equilibria of the type H + +S - ⇆ HS are must faster that other reactive processes. This means that in the stationary state of the CSTR, these reactions are actually at equilibrium, whereas these other processes only realize a collective balance. This is actually quite reasonable in the present case but could be questioned for other chemical systems. A more serious point is raised by the expression for Π ion inside the gel given by equation ( 16) since the constant K ion has to be set by reference to the values at the boundary. In a 1-D system with a CSTR at stationary state, the boundary condition is really invariant during the whole computation. Moreover, the value of this constant was found to present only small variations when the composition of the CSTR was changed within the flow state domain. Nevertheless, extension of the algorithm to non-stationary boundary conditions would need to reconsider this question. Finally, our approach is limited to the case of continuous size changes, which excludes a discontinuous volume transition. The dynamics of such transitions has only been considered on simple models with no charge, no ionic pressure and no chemical reaction [33]. At the present time, extension to such a situation is still out of reach, since, not only the dynamics of steep fronts would have to be accounted for, but the volume fraction of the polymer would be very large in the shrunken phase so that standard descriptions of transport would break down.

Fig. 2 .

 2 Fig. 2. Concentration profiles of H + and ca in the gel. F: Flow state. M: Mixed state. F/M: Spatial bistability. CSTR input flows: [BrO - 3 ]0=25 mM, [SO 2- 3 ]0=60 mM, [H2SO 2- 4 ]0=6.5 mM, l=1mm)

Fig. 3 .

 3 Fig. 3. Non equilibrium state diagram of oscillatory chemoresponsive gel. Black zone: oscillatory behavior (l is l0). Vertical bars are covering the variations of gel size l during an oscillation for three parameters sets (white points). The indicative limits of spatial bistability for fixed size are reported from Fig. 1.

Fig. 4 .

 4 Fig.4. Oscillations of size l and of [H + ] at the fixed wall for point (2) in Fig.3). [BrO - 3 ]0=25 mM, [SO 2- 3 ]0=60 mM, [H2SO4]0=11.10 mM. l0=0.25mm, D0=0.03 mm 2 /s, η=5.

Fig. 5 .

 5 Fig. 5. Oscillation shapes for points (1) and (3) in Fig. 3.

Table 2 .

 2 Diffusion coefficients. Asterisks correspond to values used for unavailable coefficients.

	Species	zi	Di × 10 -5 cm 2 s -1
	H + +1	9.312
	BrO -3 SO 2-3	-1 -2	1.485 1.1 *
	HSO -3	-1	1.5 *
	OH -	-1	5.26
	SO 2-4	-2	1.065
	HSO -4 S2O 2-6	-1 -2	1.33 1.0 *
	Br -	-1	2.084
	H2SO3	0	1.6 *
	Na + +1	1.334

Table 3 .

 3 Amplitude and period of oscillations.

	point [H2SO4]0	l0	lm	∆l	∆l/lm	T
		(mM)	(mm) (mm) (mm)	%	(h)
	(1)	11.89	0.20	0.303	0.069	23	2.54
	(2)	11.10	0.25	0.414	0.159	38	4.75
	(3)	8.65	0.385	0.730	0.313	43	17

This work has been supported by CNRS and the Agence Nationale de la Recherche. I am indebted to I. Szalai and P. De Kepper for numerous discussions on the experimental developments and communication of their most recent results.