Amélie Lelong
email: amelie.lelong@gipsa-lab.grenoble-inp.fr

Gérard Bailly
email: gerard.bailly@gipsa-lab.grenoble-inp.fr

Study of the Phenomenon of Phonetic Convergence thanks to Speech Dominoes

Keywords:

During an interaction people are known to mutually adapt. Phonetic adaptation has been studied notably for prosodic parameters such as loudness, speech rate or fundamental frequency. In most of the cases, results are contradictory and the effectiveness of phonetic convergence during an interaction remains an open issue. This paper describes an experiment based on a children game known as speech dominoes that enabled us to collect several hundreds of syllables uttered by different speakers in different conditions: alone before any interaction vs. after it, in a mediated interaction vs. in a face-to-face interaction. Speech recognition techniques were then applied to globally characterize a possible phonetic convergence.

 Introduction

The Communication Adaptation Theory (CAT), introduced by Giles et al [START_REF] Giles | Speech accommodation theory: The first decade and beyond[END_REF], postulates that individuals accommodate their communication behavior either by becoming much closer of their interlocutor (convergence) or on the contrary by increasing their differences (divergence). People can adapt to each other in different ways. For example, conversational partners notably adapt to each other"s choice of words and references [START_REF] Brennan | Lexical choice and conceptual pacts in conversation[END_REF] and also converge on certain syntactic choices [START_REF] Lockridge | Addressees needs influence speakers early syntactic choices[END_REF]. Zoltan-Ford [START_REF] Zoltan-Ford | How to get people to say and type what computers can understand[END_REF] has shown that users of dialog systems converge lexically and syntactically to the spoken responses of the system. Ward et al [START_REF] Ward | Dialog convergence and learning[END_REF] demonstrated that adaptive systems mimicking this behavior facilitate learning. This alignment [START_REF] Pickering | Activation of syntactic priming during language production[END_REF] may have several benefits such as easing comprehension [START_REF] Adank | Imitation improves language comprehension[END_REF], facilitating the exchange of messages of which the meaning is highly context-dependent [START_REF] Lakin | The chameleon effect as social glue: evidence for the evolutionary significance of nonconscious mimicry[END_REF], disclosing ability and willingness to perceive, understanding or accepting new information [START_REF] Allwood | Bodily communication -dimensions of expression and content[END_REF] and maintaining social glue or resonance [START_REF] Kopp | Social resonance and embodied coordination in face-to-face conversation with artificial interlocutors[END_REF].

Researchers have examined also adaptation of phonetic dimensions such as pitch [START_REF] Gregory | A nonverbal signal in voices of interview partners effectively predicts communication accommodation and social status perceptions[END_REF], speech rate [START_REF] Edlund | Pause and gap length in face-to-face interaction[END_REF], loudness [START_REF] Kousidis | Towards measuring continuous acoustic feature convergence in unconstrained spoken dialogues[END_REF], dispersions of vocalic targets [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] as well as more global alignment such as turn-taking [START_REF] Benus | Are we 'in sync': Turn-taking in collaborative dialogues[END_REF]. But the results of these different studies show a weak convergence and even in some cases no convergence at all. In the perceptual study conducted by Pardo [START_REF] Pardo | On phonetic convergence during conversational interaction[END_REF], disparities between talkers have been attributed to various dimensions such as social settings, communication goals and varying roles in the conversation. Sex differences have also been put forward: female interlocutors show more convergence than males.

The emerging field of research is crucial to the comprehension of adaptive behavior during unconstrained conversation on one hand and to versatile speech technologies that aim at substituting one partner with an artificial conversational agent on the other hand. Literature shows that two main challenges persist: (a) the need of original experiments that allow us to collect sufficient phonetic material to study and isolate the impact of the numerous factors influencing adaptation; (b) the use of automatic techniques for characterizing the degree of convergence if any.

 State of the art

In the following section, several influential articles will be presented. These papers thoroughly summarize research about phonetic adaptation.

Convergence and social role

There are only a few studies that explain the role of convergence in a social interaction. Different interpretations have been given.

First of all, convergence could be a consequence of the episodic memory system [START_REF] Dijksterhuis | The perception-behavior expressway: automatic effects of social perception on social behavior[END_REF]. People keep a trace of all their multimodal experiences during social interaction. An exemplar-based retrieval of previous behavior given similar social context is triggered so that the current interaction benefits from previous attunement.

Adaptation can also be used in a community to let a more stable form emerge across those present in the community [START_REF] Garrod | Conversation, co-ordination, and convention: An empirical investigation of how groups establish linguistic conventions[END_REF] or to help people to define their identity by categorizing others and themselves into groups constantly compared and evaluated [START_REF] Tajfel | An integrative theory of intergroup conflict[END_REF].

Other studies have shown that convergence may help to accomplish mutual goal [START_REF] Clark | Using Language[END_REF], align representations [START_REF] Garrod | Conversation, co-ordination, and convention: An empirical investigation of how groups establish linguistic conventions[END_REF], increase the quality of an interaction [START_REF] Babel | Phonetic and social selectivity in speech accommodation[END_REF], and furthermore contribute to mutual comprehension by decreasing social distance [START_REF] Babel | Phonetic and social selectivity in speech accommodation[END_REF].

According to Labov [START_REF] Labov | The anatomy of style-shifting, in Style and Sociolinguistic Variation[END_REF], convergence could be due to the need to add emphasis to expression and persist for the next interaction.

Finally, adaptation could be interpreted as a behavioral strategy to achieve particular social goals such as approval [START_REF] Gregory | A nonverbal signal in voices of interview partners effectively predicts communication accommodation and social status perceptions[END_REF][START_REF] Giles | Language and Social Psychology[END_REF] or desirability [START_REF] Natale | Social desirability as related to convergence of temporal speech patterns[END_REF].

Description of key studies on phonetic convergence

Pardo [START_REF] Pardo | On phonetic convergence during conversational interaction[END_REF] examined whether pairs of talkers converged in their phonetic repertoire during a single conversational interaction called a map-task. Six same-sex pairs were recruited to solve a series of 5 map tasks where their role -instruction giver or receiver -were exchanged. The advantage of the map task is to collect landmark names that are uttered several times during the interaction by each interlocutor in order to have the receiver replicate the itinerary described by the giver. One or two weeks before any interaction, talkers read out the set of map task landmark labels in order to obtain reference pronunciations. Just after interaction, the same procedure was also performed to test the persistence of convergence, i.e. to distinguish stimulusdependent mimicry from mimesis which is supposed to originate from a deeper change of phonetic representations [START_REF] Donald | Origins of the Modern Mind: three stages in the evolution of culture and cognition[END_REF]. To measure convergence, 30 listeners were asked to judge similarity between pronunciations of pre-, map-and post-task landmark labels in a AXB test, X being a map-task utterance and (A,B) pre-, map-or post-task of the same utterance pronounced by the corresponding partner. Results of this forced choice showed significant main effects of expose and persistence but there was also dependence of role and sex: givers" instructions converged more than receivers" instructions and particularly for female givers. It is in agreement with the results found by Namy [START_REF] Namy | Gender differences in vocal accommodation: The role of perception[END_REF].

Delvaux and Soquet [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] questioned the influence of ambient speech on the pronunciations of some keywords. These keywords were chosen in order to collect representatives of two sounds (the mid-open vowels [] and []) the allophonic variations of which are typical of the two dialects of French spoken in Belgium. During these non interactive experiments, subjects were asked to describe a simple scene: "C'est dans X qu'il y a N Y" (It's in X that there are N Y), where X were locations, N numbers and Y objects. This description was either uttered by the speaker or by recorded speakers using the same or the other dialect. Pre-and posttasks were also performed for the same reasons enounced previously. The phonetic analysis focused on the production of the two sounds that were used in the two possible labels X. The authors sought for unintentional imitation. To characterize the amplitude of that change, they compared durations and spectral characteristics of target sounds. In most cases, small but significant displacements towards the prototypes of the other ambient dialect were observed for both sounds (see the lowering of the canonical values in Test 1 and 2 in Fig. 1). Similar unconscious imitation of characteristics of ambient speech has also been observed by Gentilucci et al [START_REF] Gentilucci | Imitation during phoneme production[END_REF] for audiovisual stimulations. Aubanel and Nguyen [START_REF] Aubanel | Automatic recognition of regional phonological variation in conversational interaction[END_REF] also conducted experiments to study the mutual influence between French accents, i.e. northern versus southern, that could be part of the subjects" experience. They have proposed an original paradigm to collect dense interactive corpora made up of uncommon proper nouns. They defined some criteria in order to discriminate the two accents, i.e. schwa, back mid vowels, mid vowels in word-final syllables, coronal stops, and nasal vowels. Uncommon proper nouns containing these segments are chosen so as to maximize coverage of alternative spellings. They chose their subjects in a major high school and grouped them according to their sex and to a similar score on the Crowne-Marlowe [START_REF] Crowne | A new scale of social desirability independent of psychopathology[END_REF] social desirability scale. One week before any interaction, subjects read out three sets of 16 names to get reference pronunciations. This session was repeated just after the interactions to measure mimesis. During the interaction, dyads were asked to associate names with photographs and the corresponding characters" statements. Aubanel and Nguyen used a Bayes classifier to automatically assign subjects to a group and test different levels of convergence in the dyads (towards the interlocutor, the interlocutor"s group and accent) using linear discriminant analysis performed on spectral targets. They found very few instances of convergence. Additionally convergence was quite dependent of the critical segments analyzed, the sessions and the pairs.

Comments

These studies show that phonological and phonetic convergence is very weak. The experimental paradigms used so far either collect few instances (typically a dozen in Aubanel and Nguyen) of few key segments or many instances of a very small set of key segments (two in Delvaux and Soquet). These segments are always produced in a controlled context within key words.

Both studies have focused on inter-dialectal convergence and segments that carry most of the dialectal variation. This a priori choice is questionable since it remains to be shown that subjects at first negotiate these critical segments before or more easily than others. Since the convergence is segment-dependent, it is interesting to study the speakers" alignment on the common repertoire of their mother tongue. In our experiments, we will examine the convergence of the eight French peripheral oral vowels.

In most studies, interlocutors or ambient speech are not known a priori by the subjects. The authors were certainly expecting to observe on-line convergence as the dialog proceeds. The hypothesis that adaptation and alignment is immediate and fast is questionable: in the following we will compare convergence of unknowns with those of good friends.

 Material and Protocol

During our experiments, speakers were instructed to choose between two words displayed on a computer screen.

Speech Dominoes

The rule of the game is quite simple. Speakers have to choose between two words the one that begins with the same syllable as the final syllable of the word previously uttered by the interlocutor (see Table 1). Such rhyme games -here speech dominoes -are part of the children"s folklore and widely used in primary school, for example for language learning. We decided to chain simple disyllabic words such as: bateau [], taudis [], diffus [], furie [], etc. We used only two disyllabic words in order to limit the cognitive load and ease the running of successive sessions.

The words were chosen to uniformly collect allophonic variations of the eight peripheral oral vowels of French:

[], [], [], [], [], [], [], [].
To force mutual attention during the interaction, the word list has been built so that the speaker could not guess the next domino given the sole history of the dialog. In fact, he has to pay attention to the word uttered by his interlocutor to decide which "domino" he would have to utter next. For instance, spk 2, after having chosen [] in Table 1, will be presented with the following two alternatives, namely [] and []. Since [] and [] are two valid French common words with almost the same word frequency, spk 2 will have to wait until spk1 chooses the right rhyme to decide about his own.

A chain of 350 dominoes was thus established that permitted us to collect almost 40 exemplars of each peripheral oral vowel (see Table 2).

Conditions

The speakers pronounced dominoes under different conditions. First of all, we needed to get references for each speaker, we called this condition pre-test. To do this, they uttered a list of 350 words before any dialog with their interlocutor. The pre-test words were the same as those pronounced by the two speakers during the dominoes' game. It allowed us to characterize each speaker's phonetic space and to measure the amplitude of adaptation if any.

Fig. 2. Face-to-face interaction

Experiments

In this paper, we only contrast the pre-test condition and the interactive game played during three experiments:  Experiment I: speakers were in two different rooms and communicated through microphones and headphones. This setup was easy to realize thanks to the MICAL platform of our laboratory (two rooms separated with a tinted mirror). Speakers were unknown to each other.  Experiment II: same as Experiment I but with a reduced set of good friends, people that know each other or work together since a long time (mean of 15 years from 10 years to 25 years).  Experiment III: speakers were in a face-to-face interaction. We studied also dyads of good friends (mean of from 6 months to 3 year and 6 months). In the two cases, they were instructed to avoid speech overlaps and repairs so as to ease automatic segmentation and alignment.

Experimental settings

For Experiments I and II, people played through sets of microphones and headphones. Signals were digitized at 16 kHz thanks to a high-quality stereo sound card. Dominoes were displayed on two computer screens displaying a pdf file. For Experiment III, the setting is quite different. Speakers sat on each side of a table facing two back-to-back computer screens. They were recorded with a camera -a mirror allowed us to capture both interactants (see Figure 2) -their head movements were monitored using four infrared cameras (Qualisys system®).

We used two keyboards connected to the same computer to forward turns: when a speaker finishes uttering his domino, he presses a key on his keyboard to display the two choices for his next turn on his own screen.

Characterization

Delvaux and Soquet [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] noticed that global automatic analysis of spectral distributions by MFCC (Mel Frequency Cepstral Coefficients) lead to quasi-identical but more robust characterization of convergence than a more detailed semi-automatic phonetic analysis such as formant tracking. Aubanel and Nguyen [START_REF] Aubanel | Automatic recognition of regional phonological variation in conversational interaction[END_REF] similarly used automatic recognition techniques to recognize idiolects.

Here, we trained phone-sized context-independent HMMs with 5 states trained using HTK on the pre-test data. The input parameters are the first 12 MFCC + energy + deltas of theses parameters. After various forced alignments, we compared the distributions of normalized self vs. other"s recognition scores of central states of each vowel (see Fig. 4 and Fig. 5). Paired t-tests were also performed to compare changes of distributions of scores of vowels produced in the same words (175 words for each speaker).

We align signals with a network of pronunciation variants for each word to semiautomatically segment the signals. This segmentation was then checked by hand (by three different annotators). Devoiced or creaky vowels -often high vowels [], [] in unvoiced context -were discarded. Dialectal variations were also considered: allophonic variations of mid-vowels (open or closed) are determined according to the speaker-dependent partition of the range of the first formant.

 Results

Phonological variations

Despite the fact that our corpus was not designed to enhance dialectal variations (unlike Aubanel and Nguyen [START_REF] Aubanel | Automatic recognition of regional phonological variation in conversational interaction[END_REF]), we observed some dialectal variations mainly concerned with allophonic variations of mid-vowels. Most participants came from North of France and thus used exclusively open vowels in closed syllables (e.g. sabord // vs. sabot //). Other interlocutors spectrally contrasted minimal pairs such as vallée vs. valais (// vs. //), miné vs. minet (// vs. //), etc. We observed few cases of phonological adaptation, i.e. subjects adopting a pronunciation different from the one chosen in their pretest to get closer to the pronunciation of their interlocutor. Most interactions resulted in convergence of allophonic choices (see Fig. 3) but this is not significant due to limited data. For example, for the vowel [], a mutual adaptation can be seen during the interaction between ALa and MGB and also between ALa and MSM.

We should however mention that the labeling of allophonic variations of midvowels is very difficult, since French speakers have now the tendency to front midclosed vowels [START_REF] Coveney | The Sounds of Contemporary French : Articulation and Diversity[END_REF][START_REF] Boula De Mareüil | Accents étrangers et régionaux en français : Caractérisation et identification[END_REF]. The labeling is particularly difficult in non accented positions where vowel undershoot or co articulation may override perceptual intuition. We always privileged labeling based on objective measurements that tend to favor midclosed options.

Fig. 3. Proportion of peripheral mid vowels pronounced as closed by 4 pairs (left: the vowel [e]; right: vowel [o]). The initiator was the same female ALa interacting with 3 males (MGB, MMP, MSM) and one female (FLD). For each pair, bars represent the proportion uttered during respectively the ALa pretest, ALa interacting with her interlocutor, her interlocutor interacting with ALa and the interlocutor's pretest.

Sub-phonemic convergence

Given the assumption that each phone was properly labeled, we compared the distributions of normalized recognition scores of the pre-test and interactive utterances, as explained previously. These utterances were recognized by the speaker"s own HMMs in a first time and in a second time by the HMMs of his interlocutor. For pre-test data, we expected high scores for HMMs tested on their own training data by construction and lower scores for HMMs of the interlocutor. The recognition score of each vowel is the average log likelihood per frame for the central state of the corresponding HMM. The difference between the scores somehow reflects the inter-speaker distance. Convergence would be characterized by a decrease of scores by self HMMs and an increase of scores by the other's HMMs.

The recognition is thus performed by HMM models of each speaker and of his/her interlocutor. Fig. 4 and Fig. 5 compare the distributions of normalized recognition scores for the pre-test (left) versus the interaction (right). Scores are typically higher for phones uttered by one speaker and recognized by his own HMMs. In case of an interaction between two unknowns (cf. Fig. 4), the distributions computed for the interactive speech do not evolve so much. We observe stronger convergence in case of good friends (cf. Fig. 5).

Fig. 4. Distribution of recognition scores for the vowels of disyllabic words produced by two unknowns.

The recognition is performed by their own HMM models and by the HMM models of their interlocutor. Scores are expected to be higher by using their own HMM models. Left: scores for word lists read aloud in isolation; this speech data is used to train the speaker-specific HMM models. It can be seen that, by using the own model of each interlocutor (lmb_lmb and rl_rl), higher recognition scores are obtained than for cross recognition (lmb_rl and rl_lmb) Right: same words pronounced in a verbal domino game. In this case, we expect a decrease of recognition scores by using the own model of each interlocutor (lmb_rl_lmb and rl_lmb_rl) and an increase of recognition scores by using cross recognition (lmb_rl_rl and rl_lmb_lmb). Here, only small adjustments are observable (weak shift on the left for lmb_rl_lmb and rl_lmb_rl and on the right for lmb_rl_rl and rl_lmb_lmb).

Distributions of recognition scores

Fig. 6 shows the average convergence rate for all dyads recorded in the two experiments. This is computed as the relative distance of vocalic targets produced in pretest vs. interaction (central state of the HMM alignment). Similarly to Delvaux & Soquet [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF], linear discriminant analysis is performed on the target MFCC parameters for each vowel to categorize each interlocutor's vocalic space into two distinct groups. For each pair, pre-test and interactive vocalic targets are projected onto the first discriminant axis. A normalized convergence rate is then computed by dividing the distance between targets produced during the interaction with that produced during the pre-test for the same word. Convergence is not systematic. We can see that the phenomenon is amplified with pairs of the same sex and particularly women (it can be due to a priori more similar speakers). We used this statement for the seven last interactions by selecting only women for the experiment and the results confirm this statement. An ANOVA has been done to assess significance of adaptation. Distributions with significant convergence rates are drawn in bold.

Convergence of vocalic targets

Convergence is a vowel-and interlocutor-dependent phenomenon. Fig. 7 shows that some pairs do not adapt at all while others show a significant mutual adaptation.

To study interlocutor-specific adaptation strategies, our game initiators interacted with 2 to 5 different interlocutors. Fig. 8 illustrates examples of target-and interlocutor-specific behaviors. When considering each vowel separately (40 occurrences in average, see Table 1), we do observe cases of full convergence (speakers A and B in Fig. 8).

Note however that our analysis is based on relative distances and should take into account the whole structure of the vocalic space of the interlocutors. Speakers notably fill differently their available acoustic space, especially between mid-vowels [START_REF] Ménard | Invariance and variability in the production of the height feature in French vowels[END_REF].

An evolution of convergence rates with time was expected. Convergence rates as a function of time as been plotted for each interlocutor but nothing relevant has been observed. Maybe the proposed task was too short to observe this phenomenon.

Prosody

Fig. 10 shows that fundamental frequency register was relatively unaffected by the interaction. The exchange of simple words does not favor attunement of melody.

Convergence of speech rhythm is clearly observed, certainly due to the "ping-pong" task. This can also be due to the fact that speech rhythm was also much quicker in interactive speech compared to isolated word reading (cf. Fig. 10) with a notable shortening of final syllables. This is probably due to the task focusing on rhyme matching. Delvaux and Soquet [START_REF] Delvaux | The influence of ambient speech on adult speech productions through unintentional imitation[END_REF] advise to discard final syllables for studying phonetic convergence. In our case, we did not find any difference between global and partial statistics except stronger convergence for the durations of vocalic nuclei of initial syllables.

 Conclusions and perspectives

We proposed here an original speech game that quickly collects many instances of target sounds with a mutual influence that force interlocutors to engage into active action-perception loops. Distribution of target sounds can be explicitly controlled to observe convergence in action if any.

We found occurrences of strong phonetic convergence with only one instance of small divergence. This convergence strongly depended on the dyads -with strongest convergence observed for pairs of the same sex -and seemed to be phonemedependent. We used this observation to select our last subjects and the results confirmed a strongest convergence for dyads composed of women.

These objective measurements should be confirmed by subjective assessments such as promoted by Pardo [START_REF] Pardo | On phonetic convergence during conversational interaction[END_REF]. We are also planning to conduct a series of subjective tests to determine if adapted stimuli offer a clearer perceptual benefit for listeners compared as to non adapted stimuli. Perception of degraded stimuli such as used by Adank et al [START_REF] Adank | Imitation improves language comprehension[END_REF] is an interesting option.

This gaming paradigm will now be used to select subjects and dyads who exhibit the strongest adaptation abilities and study more complex conversational situations. This data will be used to train speech synthesis engines that will implement these adaptation strategies. Such interlocutor-aware components are certainly crucial for creating social rapport between humans and virtual conversational agents [START_REF] Gratch | Creating rapport with virtual agents[END_REF].

Fig. 1 .

 1 Fig. 1. Results on spectral distance calculated by Delvaux and Soquet [14]. It can be seen that, during tests (Tests 1 & 2), subjects are getting away from their own reference (Pretest) and closer to the other dialect (References 1 & 2).

Fig. 5 .

 5 Fig. 5. Same as Fig. 4 but for disyllabic words produced by two old friends. Stronger convergence is observed here since larger shifts are observed.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Average convergence rate (calculated on 100 iterations) of vocalic targets of interlocutors for all conditions. A linear discriminant analysis on one random half of the pre-test data has been used to separate each interlocutor's vocalic space. Thus, a reference discriminant distance between interlocutors is obtained. It is used to calculate normalized convergence rates. First, it is used to calculate the convergence rate on the other half of the pre-test to have our reference departure for each interlocutor. The two dotted lines represent the mean of pre-test of the tested subjects (line 0) and of the reference subjects (line 1). Distributions displayed with bold lines are significantly different (p<0.05) from the corresponding pre-test (reference departure). Note that only two significant divergences are found (one speaker in the pair number 12 and one on the pair 22). Most convergence cases are observed with pairs with same sex.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Mean changes of relative difference between F0 registers. As expected the values are higher for pairs with different sex. No significant narrowing of this difference is induced by interaction. This can be due to the task that imposes short utterances.

Table 1 . First speech dominoes used in the interactive scenario. Interlocutors have to choose and utter alternatively the rhyming words. Correct chainings of rhymes are high lightened with a dark background.

 1

	spk 1	spk 2	spk 1	spk 2	spk 1	spk 2	spk 1
	 	 	 	 	 	 	 	…

Table 2 . Number of phones collected for each speaker during the dominoes' game. 350 CV or CVC syllables are pronounced in total.

 2

	phones								 others
	#items	47	48	45	43	44	40	43	31	9

Acknowledgments

This work has been financed by ANR Amorces and by the Cluster RA ISLE. We thank Frederic Elisei, Sascha Fagel and Loïc Martin for their help.


N Mellouli
email: mellouli@math.univ-lyon1.fr

Projectively Equivariant Quantization and Symbol calculus in dimension 1|2

The spaces of higher-order differential operators (in Dimension 1|2), which are modules over the stringy Lie superalgebra K(2), are isomorphic to the corresponding spaces of symbols as orthosymplectic modules in non resonant cases. Such an osp (2|2)-equivariant quantization, which has been given in second-order differential operators case, keeps existing and unique. We calculate its explicit formula that provides extension in particular order cases.

Introduction and the main results

Let S 1|2 be a supermanifold which is endowed with a projective structure (Susy-structure), see [START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF], in dimension 1|2 and D λ,µ S 1|2 the space of differential operators on S 1|2 acting from the space of λ-densities to the space of µ-densities where λ and µ are real or complex numbers. The space D λ,µ S 1|2 which is a module over the stringy superalgebra K(2), see [START_REF] Grozman | Lie superalgebras of string theories[END_REF], is naturally filtrated and has an other finer filtration by the contact order of differential operators. The space of symbols S S 1|2 , that is the graded module grD λµ S 1|2 , isn't isomorphic to the space D λ,µ S 1|2 as K(2)-module. Therefore, we have restricted the module structure on D λµ S 1|2 to the orthosymplectic Lie superalgebra osp (2|2) that is naturally embedded into K [START_REF] Conley | Conformal symbols and the action of contact vector fields over the superline[END_REF]. We establish a canonical isomorphism between the space of differential operators on S 1|2 and the corresponding space of symbols. An explicit expression of projectively equivariant quantization map is given in case of second order differential operetors, see [START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF]. We extand calculus to symbols of higher order differential operators.

Geometry of the supercircle S 1|2

We have considered the supercircle S 1|2 described in [START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF] by its graded commutative algebra of complex-valued functions C ∞ S 1|2 , consisting of the following elements :

f (x, ξ 1 , ξ 2 ) = f 0 (x) + ξ 1 f 1 (x) + ξ 2 f 2 (x) + ξ 1 ξ 2 f 12 (x), (2.1)
where x is the Fourier image of the angle parameter on S 1 , ξ 1 , ξ 2 are odd Grassmann coordinates and f 0 , f 12 , f 1 , f 2 ∈ C ∞ (S) are functions with complex values. We have defined the parity function p by setting p (x) = 0 and p (ξ i ) = 1. The standard contact structure on S 1|2 , known as Susy-structure, is defined by the data of a linear distribution D 1 , D 2 on S 1|2 generated by the odd vector fields :

D 1 = ∂ ξ 1 -ξ 1 ∂ x , D 2 = ∂ ξ 2 -ξ 2 ∂ x . (2.2)
We would rather recall that every contact vector field can be expressed, for some function f ∈ C ∞ S 1|2 , by

X f = f ∂ x -(-1) p(f ) 1 2 D 1 (f ) D 1 + D 2 (f ) D 2 (2.
3)

The projective (conformal) structure on the supercircle S 1|2 , see [START_REF] Ovsienko | Classification of contact-projective structures on the supercircle[END_REF], is defined by the action of Lie superalgebra osp (2|2). The orthosymplectic algebra osp (2|2) is spanned by the contact vector fields X f with the contact Hamiltonians f which are elements of

{1, ξ 1 , ξ 2 , x, ξ 1 ξ 2 , xξ 1 , xξ 2 , x 2 }. The embedding of osp (2|2) into K (2) is given by (2.3).
The subalgebra Aff (2|2) of osp (2|2), called the Affine Lie superalgebra, is spanned by the contact vector fields X f with the contact Hamiltonians f which are elements of {1, ξ 1 , ξ 2 , x, ξ 1 ξ 2 }.

For any contact vector field, we have defined a family of differential operators of order one on

C ∞ S 1|2 : L λ X f := X f + λf ′ , (2.4)
where the parameter λ is an arbitrary complex number. Thus, we have obtained a family of K (2)-modules on C ∞ S 1|2 noted by F λ S 1|2 which are called the spaces of weighted densities of weight λ.

Differential operators on the spaces of weighted densities

In this section, we have introduced the space of differential operators acting on the spaces of weighted densities. We have also presented the corresponding space of symbols on S 1|2 . Those are detailed in [START_REF] Leites | New invariant differential operators on supermanifolds and pseudo-(co)homology[END_REF][START_REF] Grozman | Invariant operators on supermanifolds and standard models[END_REF][START_REF] Gargoubi | Differential operators on supercircle: conformally equivariant quantization and symbol calculus[END_REF][START_REF] Conley | Conformal symbols and the action of contact vector fields over the superline[END_REF][START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF]. For every integer or half-integer k, the space of differential operators of the form

A = ℓ+ m 2 + n 2 ≤k a ℓ,m,n ∂ ℓ x D m 1 D n 2 , (3.1)
where a ℓ,m,n ∈ C ∞ S 1|2 and m, n ≤ 1, has been noted by D k λµ S 1|2 . This above K(2)module space has a K(2)-invariant finer filtration :

D 0 λµ S 1|2 ⊂ D 1 2 λµ S 1|2 ⊂ D 1 λµ S 1|2 ⊂ ... ⊂ D k λµ S 1|2 ⊂ D k+ 1 2 λµ S 1|2 ⊂ • • • (3.2)
that has been considered in papers [START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF][START_REF] Gargoubi | Differential operators on supercircle: conformally equivariant quantization and symbol calculus[END_REF].

We would rather remind that the orthosymplectic superalgebra osp (2|2) has been acting on D k λµ S 1|2 . The action of contact field X f of order one on a differential operator A, see formula (2.3), is given by

£ λµ X f (A) = L µ X f • A -(-1) p(f )p(A) A • L λ X f . (3.3)

Space of symbols of differential operators

The graded K (2)-module, which is associated with the finer filtration (3.2) and called the space of symbols of differential operators, is defined by

grD λµ S 1|2 = ∞ i=0 gr i 2 D λµ S 1|2 , (3.4)
where

gr k D λµ S 1|2 = D k λµ S 1|2 /D k-1 2 λµ S 1|2
for every integer or half-integer k. The image of any differential operator through the natural projection

σ pr : D k λµ S 1|2 → gr k D λµ S 1|2
, that is defined by the filtration (3.2), has been called the principal symbol.

Referring to [START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF], the stringy superalgebra K (2) keeps acting on the space of symbols.

Proposition 3.1. If k is an integer, then

gr k D λµ S 1|2 = F µ-λ-k F µ-λ-k (3.5)
Proof. By definition (see formula (3.1)), a given operator

A ∈ D k λµ S 1|2 with integer k is of the form A = F 1 ∂ k x + F 2 ∂ k-1 x D 1 D 2 + • • •
, where • • • stand for lower order terms. The principal symbol of A is then encoded by the pair (F 1 , F 2 ). From (3.3), we easily calculate the K(2)-action on the principal symbol:

L X f (F 1 , F 2 ) = L µ-λ-k X f (F 1 ) , L µ-λ-k X f (F 2 ) .
In other words, both F 1 and F 2 transform as (µ -λ -k)-densities.

The situation is more complicated for half-integer k : the K (2)-action has been given by

L X f (F 1 , F 2 ) = (L µ-λ-k X f (F 1 ) - 1 2 D 1 D 2 (f ) F 2 , L µ-λ-k X f (F 2 ) + 1 2 D 1 D 2 (f ) F 1 ). (3.6)
Therefore, the spaces of symbols of half-integer contact order aren't isomorphic to the spaces of weighted densities. Simplifying the notation as in [START_REF] Lecomte | Projectively equivariant symbol calculus[END_REF][START_REF] Duval | Conformally equivariant quantization: existence and uniqueness[END_REF][START_REF] Gargoubi | Differential operators on supercircle: conformally equivariant quantization and symbol calculus[END_REF][START_REF] Mellouli | Second-order Conformally Equivariant Quantization in Dimension 12[END_REF], the whole space of symbols grD λµ S 1|2 , depending only on µ -λ, has been noted by S µ-λ S 1|2 , and the space of symbols of contact order k has been noted by S k µ-λ S 1|2 . A linear map, Q : S µ-λ S 1|2 → D λµ S 1|2 , is called a quantization map if it verifies bijectivety and preserves the principal symbol of every differential operator.

Projectively equivariant quantization on S 1|2

The main result of this paper is the existence and uniqueness of an osp (2|2)-equivariant quantization map in Dimension 1|2. We calculate its explicit formula.

For every m integer or half-integer, the space D m λµ S 1|2 is isomorphic to the corresponding space of symbols as an Aff (2|2)-module. We will show how to extend this isomorphism to that of the osp (2|2)-modules.

The Divergence operators as Affine equivariant

Let us introduce new differential operators, which are called Divergence operators, on the space of symbols S µ-λ S 1|2 .

At first, we consider the case of differential operators of contact order k, where k is an integer. We have assumed that the symbols of differential operators are homogeneous and we have defined parity of the non vanished symbol (F 1 , F 2 ) as p (F ) := p (F 1 ) = p (F 2 ).

The Divergence operators in case of integer contact order k

In this case, we define the Divergence as Affine equivariant differential operators on the space of symbols S µ-λ S 1|2 . In each component S k µ-λ S 1|2 , we have

DIV 2n+1 (F 1 , F 2 ) = (-1) p(F )+1 k+2λ k ∂ n x D 2 (F 2 ) + ∂ n x D 1 (F 1 ) ∂ n x D 2 (F 1 ) -k+2λ k ∂ n x D 1 (F 2 ) , (4.1)
DIV 2n (F 1 , F 2 ) = ∂ n x (F 1 ) - (k+2λ)n k(2(µ-λ)+n-2k) ∂ n-1 x D 1 D 2 (F 2 ) (k+2λ)(k-n) k(k-n+2λ) ∂ n x (F 2 ) + n(k-n) (2(µ-λ)+n-2k)(k-n+2λ) ∂ n-1 x D 1 D 2 (F 1 ) (4.2)
and

div 2k-(2n+1) = ∂ k-n-1 x D 1 , ∂ k-n-1 x D 2 , div 2k-(2n) = ∂ k-n x , ∂ k-n-1 x D 1 D 2 .
Lemma 4.1. The Divergence operators (4.1) and (4.2) commute with the Aff (2|2)-action.

Proof. This is a direct consequence of projectively equivariant symbol calculus. We are looking for the symbols:

DIV 2n+1 (F 1 , F 2 ) = c 1 ∂ n x D 2 (F 2 ) + c 2 ∂ n x D 1 (F 1 ) c 3 ∂ n x D 2 (F 1 ) + c 4 ∂ n x D 1 (F 2 )
and

DIV 2n (F 1 , F 2 ) = c 5 ∂ n x (F 1 ) + c 6 ∂ n-1 x D 1 D 2 (F 2 ) c 7 ∂ n x (F 2 ) + c 8 ∂ n-1 x D 1 D 2 (F 1 )
where c i (1 ≤ i ≤ 8) are arbitrary constants. From the commutation relation L X f , DIV for f ∈ Aff (2|2) we easily get the Aff (2|2)-equivariance if Divergence operators . In this case, we also define the Divergence as Affine equivariant differential operators on the space of symbols S µ-λ S 1|2 . In each component S

k+ 1 2 µ-λ S 1|2 we have DIV 2n+1 (F 1 , F 2 ) = (-1) p(F ) 2 (µ -λ) -(2k + 1) 2 (µ -λ) -2k ∂ n x D 2 (F 2 ) + ∂ n x D 1 (F 1 ) k-n k-n+2λ ∂ n x D 2 (F 1 ) -∂ n x D 1 (F 2 ) , (4.3)
DIV 2n (F 1 , F 2 ) = 2 (µ -λ) -(2k + 1) 2 (µ -λ) -2k 2(µ-λ)+n-2k 2(µ-λ)+n-(2k+1) ∂ n x (F 1 ) - n 2(µ-λ)+n-(2k+1) ∂ n-1 x D 1 D 2 (F 2 ) 2(µ-λ)+n-2k 2(µ-λ)+n-(2k+1) ∂ n x (F 2 ) + n 2(µ-λ)+n-(2k+1) ∂ n-1 x D 1 D 2 (F 1 ) , (4.4
) and

div 2k+1-(2n+1) = ∂ k-n x , ∂ k-n-1 x D 1 D 2 , div 2k+1-(2n) = ∂ k-n x D 1 , ∂ k-n x D 2 .
Lemma 4.2. The Divergence operators (4.3) and (4.4) commute with the action of Affine Lie superalgebra.

Proof. Straightforward calculus.

Staitement of the main result

Let us give the explicit formula of the projectively equivariant quantization map. We will give the proof in the next section.

Theorem 4.3. The unique osp (2|2)-equivariant quantization map associates the following differential operator with a symbol

(F 1 , F 2 ) ∈ S k µ-λ S 1|2
where k is (even or odd) integer :

Q (F 1 , F 2 ) = k n=0 k 2 2n+1+(-1) k 4 k-1 2 + 2λ 2n+1-(-1) k 4 k -2 (µ -λ) n+1 2 DIV n (F 1 , F 2 ) div k-n (4.5) provided µ -λ = 0, 1 2 , 1, 3 2 , 2.
..where DIV and div are defined in each particular case of even or odd contact order : (4.1), (4.2), (4.3) and (4.4) and the coefficients are

n m = n(n-1)..(n-m+1) m! .
Remark 4.4. This theorem keeps being achieved in the case of Dimension 1|1, see [START_REF] Gargoubi | Differential operators on supercircle: conformally equivariant quantization and symbol calculus[END_REF]; the divergence operators DIV and div are given by D.

Proof of theorem in case of k-differential operators

Proof. Let us first consider the case of k-differential operators where k is integer. The quantization map 4.5 is, indeed, osp (2|2)-equivariant. Now, we are considering a differentiable linear map Q : S k µ-λ S 1|2 → D k λµ S 1|2 for k ≥ 1, preserving the principal symbol. Such a map is of the form :

Q (F 1 , F 2 ) = F 1 ∂ k x + F 2 ∂ k-1 x D 1 D 2 + ... + Q (2ℓ) 1 (F 1 ) + Q (2ℓ+1) 1 (F 1 ) + Q (2ℓ) 2 (F 2 ) + Q (2ℓ+1) 2 (F 2 ) .. + C 2k,1 ∂ k x (F 1 ) + C 2k,2 ∂ k-1 x D 1 D 2 (F 2 ) where Q (m) 1 and Q (m) 2
are differential operators with coefficients in F µ-λ S 1|2 , see (3.1). We obtain the following : a) This map commutes with the action of the vector fields D 1 , D 2 ∈ osp (2|2), where

D i = ∂ ξ i + ξ i ∂ x , if and only if the differential operators Q (m) 1 and Q (m) 2
are with constant coefficients.

b) This map commutes with the linear vector fields X ξ 1 , X ξ 2 , X x if and only if the differential operators

Q (m) 1 and Q (m) 2 are of contact order m 2 in addition to the form              Q (2ℓ+1) 1 (F 1 ) = C 2ℓ+1,1 ∂ ℓ x D 1 (F 1 ) ∂ k-ℓ-1 x D 1 + C 2ℓ+1,3 ∂ ℓ x D 2 (F 1 ) ∂ k-ℓ-1 x D 2 Q (2ℓ) 1 (F 1 ) = C 2ℓ,1 ∂ ℓ x (F 1 ) ∂ k-ℓ x + C 2ℓ,3 ∂ ℓ-1 x D 1 D 2 (F 1 ) ∂ k-ℓ-1 x D 1 D 2 Q (2ℓ+1) 2 (F 2 ) = C 2ℓ+1,2 ∂ ℓ x D 2 (F 2 ) ∂ k-ℓ-1 x D 1 + C 2ℓ+1,4 ∂ ℓ x D 1 (F 2 ) ∂ k-ℓ-1 x D 2 Q (2ℓ) 2 (F 2 ) = C 2ℓ,2 ∂ ℓ-1 x D 1 D 2 (F 2 ) ∂ k-ℓ x + C 2ℓ,4 ∂ ℓ x (F 2 ) ∂ k-ℓ-1 x D 1 D 2
where the coefficients C m,i (i = 1, 2, 3, 4) are arbitrary constants. Note that the vector field X xξ i is the commutation relation [X ξ i , X x 2 ] , i = 1, 2, so it is sufficient to impose the equivariance with respect to the vector field X x 2 to meet the whole condition of osp (2|2)-equivariance.

d) The above quantization map commutes with the action of X x 2 if and only if any the coefficients C m,i verify the following conditions :

                                           ℓ (ℓ -1 + 2 (µ -λ -k)) C 2ℓ,1 = -(k -ℓ + 1) (k -ℓ + 2λ) C 2ℓ-2,1 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ,2 = (k -ℓ + 2λ) (-1) p(F ) (C 2ℓ-1,2 -C 2ℓ-1,4 ) -(k -ℓ + 1) C 2ℓ-2,2 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ,3 = -(k -ℓ) (-1) p(F ) (C 2ℓ-1,1 + C 2ℓ-1,3 ) + (k -ℓ + 2λ + 1) C 2ℓ-2,3 ℓ (ℓ -1 + 2 (µ -λ -k)) C 2ℓ,4 = -(k -ℓ) (k -ℓ + 2λ + 1) C 2ℓ-2,4 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ+1,1 = (k -ℓ) (-1) p(F ) C 2ℓ,1 -(k -ℓ + 2λ) C 2ℓ-1,1 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ+1,2 = (k -ℓ + 2λ) (-1) p(F ) C 2ℓ,4 -(k -ℓ) C 2ℓ-1,2 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ+1,3 = (k -ℓ) (-1) p(F ) C 2ℓ,1 -(k -ℓ + 2λ) C 2ℓ-1,3 (ℓ + 1) (ℓ + 2 (µ -λ -k)) C 2ℓ+1,4 = -(k -ℓ + 2λ) (-1) p(F ) C 2ℓ,4 + (k -ℓ) C 2ℓ-1,4
If µ-λ = 0, 1 2 , 1, 3 2 , 2..., this system has been solved and the solutions are the following:

                       C 2ℓ,2 =   k -1 ℓ -1     k + 2λ ℓ + 1     -2 (µ -λ -k) ℓ + 1   C 2ℓ,3 = -   k ℓ + 1     k + 2λ -1 ℓ -1     -2 (µ -λ -k) ℓ + 1   and                    C 2ℓ,1 = (k+2λ-ℓ)(2(µ-λ-k)+ℓ) ℓ(k-ℓ) C 2ℓ,3 , C 2ℓ,4 = -(k-ℓ)(2(µ-λ-k)+ℓ) ℓ(k-ℓ+2λ) C 2l,2 , C 2ℓ+1,1 = (-1) p(F ) (k+2λ-ℓ) ℓ C 2ℓ,3 , C 2ℓ+1,2 = -(-1) p(F ) (k-ℓ) ℓ C 2ℓ,2 , C 2ℓ+1,3 = (-1) p(F ) (k+2λ-ℓ) ℓ C 2ℓ,3 C 2ℓ+1,4 = (-1) p(F ) (k-ℓ) ℓ C 2ℓ,2
.

That allows us to obtain the formula (4.5).

Proof of theorem in case of k + 1 2 -differential operators

Proof. In the case of k + 1 2 -differential operators where k is integer, we get an Aff (2|2)equivariant quantization map by a straightforward calculation which is given by

Q (F 1 , F 2 ) = F 1 ∂ k x D 1 + F 2 ∂ k x D 2 + ... + Q (2ℓ) 1 (F 1 ) + Q (2ℓ+1) 1 (F 1 ) + Q (2ℓ) 2 (F 2 ) + Q (2ℓ+1) 2 (F 2 ) .. + C 2k+1,1 ∂ k x D 1 (F 1 ) + C 2k+1,2 ∂ k x D 2 (F 2 )
where the m 2 -differential operators

Q (m) 1 and Q (m) 2
have the form :

             Q (2ℓ) 1 (F 1 ) = C 2ℓ,1 ∂ ℓ x (F 1 ) ∂ k-ℓ x D 1 + C 2ℓ,3 ∂ ℓ-1 x D 1 D 2 (F 1 ) ∂ k-ℓ x D 2 Q (2ℓ+1) 1 (F 1 ) = C 2ℓ+1,1 ∂ ℓ x D 1 (F 1 ) ∂ k-ℓ x + C 2ℓ+1,3 ∂ ℓ x D 2 (F 1 ) ∂ k-ℓ-1 x D 1 D 2 Q (2ℓ) 2 (F 2 ) = C 2ℓ,2 ∂ ℓ-1 x D 1 D 2 (F 2 ) ∂ k-ℓ x D 1 + C 2ℓ,4 ∂ ℓ x (F 2 ) ∂ k-ℓ x D 2 Q (2ℓ+1) 2 (F 2 ) = C 2ℓ+1,2 ∂ ℓ x D 2 (F 2 ) ∂ k-ℓ x + C 2ℓ+1,4 ∂ ℓ x D 1 (F 2 ) ∂ k-ℓ-1 x D 1 D 2
The above quantization map commutes with the action of X x 2 if and only if the coefficients C m,j (j = 1, 2, 3, 4) verify the following system of linear equations : By solving this system, we obtain the formula (4.5).

                                                                                           C 2ℓ,1 -(ℓ + 1) ℓ + 2 µ -λ -k -

4. 1 . 2

 12 The Divergence operators in case of half-integer contact order k + 1 2

2 C 1 C 2 C C 2ℓ+1, 1 + 2 C 2ℓ+1, 2  2 C 2ℓ+1, 1  C 2ℓ+1, 3 --1 2 C 2ℓ+1, 4  C 2ℓ+1, 4 --1 2 C 2ℓ+1, 3 

 21212221324423 ℓ + 1) (k -ℓ + 2λ + 1) C 2ℓ-2,2 + (-1) p(F ) (k -ℓ + 1) C 2ℓ-1,2 +C 2ℓ,4 -(-1) p(F ) (k -ℓ + 2λ + 1) C 2ℓ-11) p(F ) (k -ℓ + 1) C 2ℓ-1,1 -(k -ℓ + 1) (k -ℓ + 2λ + 1) C 2ℓ-2,3 + (-1) p(F ) (k -ℓ + 2λ + 1) C 2ℓ-11 + 2 µ -λ -k -1 2ℓ,1 -C 2ℓ,2 = -(k -ℓ + 1) (k -ℓ + 2λ + 1) C 2ℓ-2,2ℓ,3 +ℓ ℓ -1 + 2 µ -λ -k -1 2ℓ,4 = -(k -ℓ + 1) (k -ℓ + 2λ + 1) C 2ℓ-2,4  (ℓ + 1) ℓ + 2 µ -λ -k -1 = -(k -ℓ + 2λ) (-1) p(F ) C 2ℓ,4 + (k -ℓ + 1) C 2ℓ-1,2   C 2ℓ+1,2 + (ℓ + 1) ℓ + 2 µ -λ -k -1 = -(k -ℓ + 2λ) (-1) p(F ) C 2ℓ,1 + (k -ℓ + 1) C 2ℓ-1,1  (ℓ + 1) ℓ + 2 µ -λ -k = (k -ℓ) (k -ℓ + 2λ + 1) C 2ℓ-1,4 -(-1) p(F ) C 2ℓ,4  (ℓ + 1) ℓ + 2 µ -λ -k = (k -ℓ) (-1) p(F ) C 2ℓ,1 + (k -ℓ + 2λ + 1) C 2ℓ-1,3

Acknowledgements. I am very grateful to C.Duval, H. Gargoubi and particularly to my supervisor V. Ovsienko for providing me with the problem and the constant help. I am also thankful to D. Leites for his critical reading of this paper and his helpful suggestions.


François Boulier

François Lemaire

Michel Petitot

Alexandre Sedoglavic

Chemical Reaction Systems, Computer Algebra and Systems Biology

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

Introduction

In this invited paper, we survey some of the results obtained in the computer algebra team of Lille, in the domain of systems biology. So far, we have mostly focused on models (systems of equations) arising from generalized chemical reaction systems. Eight years ago, our team was involved in a joint project, with physicists and biologists, on the modeling problem of the circadian clock of the green algae Ostreococcus tauri. This cooperation led us to different algorithms dedicated to the reduction problem of the deterministic models of chemical reaction systems. More recently, we have been working more tightly with another team of our lab, the BioComputing group, interested by the stochastic dynamics of chemical reaction systems. This cooperation led us to efficient algorithms for building the ODE systems which define the statistical moments associated to these dynamics. Most of these algorithms were implemented in the MAPLE computer algebra software. We have chosen to present them through the corresponding MAPLE packages.

Chemical reaction systems provide a general setting for modeling in biology. More generally, chemical kinetics may be viewed as a prototype of nonlinear science, as pointed out in [33, chapter 1].

The reaction A ----→ B describes the transformation of a species A into a different species B. Species A is the reactant of the reaction. Species B is the product. One often endows a reaction with a symbol k, which parametrizes the speed of the transformation. The reaction is then denoted A k ----→ B. The reaction A ----→ ∅ describes the transformation of species A into a species which is not part of the model. This sort of reaction often occurs when one models biological phenomenons but one usually does not encounter it in chemistry, since it is not equilibrated. Symmetrically, the reaction ∅ ----→ B describes the entry, in the model, of a species B, from outside the model. A more complicated reaction is A + B ----→ C. It is interpreted as follows: when a molecule of A encounters a molecule of B, both may react and form a molecule of a third species C. Last, one sometimes encounters reactions denoted A + B ----→ ← ----C. Some authors consider it as a single revertible reaction. We view it as a pair of two reactions. The right to the left reaction may be described as follows: every molecule of C may break itself and yield two molecules: one of species A and one of species B.

A chemical reaction system is a set of chemical reactions. Here are two examples. The first one is a classics of chemistry lectures: it is the simplest example of an enzymatic reaction. It describes the transformation of a substrate S into a product P , in the presence of some enzyme E. An intermediate complex ES is formed:

E + S k1 ----→ ← ---- k-1 ES k2 ----→ E + P . (1)
The second example, adapted from [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF]Syst. (7.19)], models a two-species oscillator:

∅ a ----→ A , A k1 ----→ ∅ , 2 A + B k2 ----→ 3 A , ∅ b ----→ B . (2)
In order to build systems of equations, i.e. mathematical models of chemical reaction systems, one needs more precise assumptions. There are at least two families of models: the deterministic and the stochastic ones.

Deterministic modeling

There are different ways to associate a precise deterministic dynamics to a chemical reaction system. In this paper, we focus on the mass-action dynamics. The parameters of the reactions are called kinetic constants. To each species A, one associates a function A(t) which represents the concentration of the species. The evolution of each concentration is given by an ordinary differential equation, in which kinetic constants appear as parameters. The system of ordinary differential equations built using the mass-action law, from a chemical reaction system, is called the natural deterministic model of the system. It is given by the following formula:

dX dt = N • V
where X is the column vector of the concentrations, N is the stoichiometry matrix, and V is the column vector of the laws. The law of a reaction is obtained by multiplying the kinetic constant by the product of the concentrations of the reactants. The stoichiometry matrix involves one row per species and one column per reaction. Its coefficient, row r and column c, is equal to the number of molecules of species r produced by the reaction c. This number is equal to the number of occurences of species r as a product, minus the number of occurences of species r as a reactant. On System (2), we have

X = A(t) B(t) , N = 1 -1 1 0 0 0 -1 1 , V =     a k 1 A(t) k 2 A(t) 2 B(t) b     .
Natural deterministic models of chemical reaction systems were intensively studied. Many informations can be read in the stoichiometry matrix alone: the nullspace of its transpose provides linear conservation laws. System (2) does not have any, while computing a basis of this nullspace for System (1) yields the two laws:

-E(t) + S(t) + P (t) = cst 1 , E(t) + ES(t) = cst 2 . (3)
A striking property of chemical reaction systems is a very simple necessary and sufficient condition for a system of polynomial ODE to be the natural deterministic model of a chemical reaction system: it is necessary and sufficient that, in the right hand side of the ODE which describes the evolution of any concentration A(t), there does not exist any monomial which is both free of A(t) and endowed with a minus sign [33, section 4.7.1]. Another striking, but much more difficult result is the zero deficiency theorem of Feinberg, Horn and Jackson [START_REF] Feinberg | The Existence and Uniqueness of Steady States for a Classe of Chemical Reaction Networks[END_REF], which gives a sufficient condition for a system to admit a unique steady state with strictly positive coordinates. This sufficient condition, which can be formulated in the setting of the graph theory, can be implemented very easily. This theorem was much studied, recently, from different points of views. See [START_REF] Chaves | State-Estimators for Chemical Reaction Networks of Feinberg-Horn-Jackson Zero Deficiency Type[END_REF][START_REF] Gatermann | Toric ideals and graph theory to analyze Hopf bifurcations in mass action systems[END_REF] and the references therein.

The size of natural deterministic models becomes so large, on all but the simplest examples, that it forbids any further analysis. The computer algebra team of Lille has focused on the problem of reducing models. We make use of two types of reduction: a reduction method which permits to approximate the natural deterministic model, under some simplifying assumptions ; and a reduction method, which permits to reduce the number of parameters, and to separate parameters which have an effect on the coordinates and the stabilities of the steady points, from the ones which have an effect on their stabilities only.

Approximating Models

The approximation method implements the well-known quasi-steady state approximation technique. See [START_REF] Kokotovic | Singular Perturbation Methods in Control: Analysis and Design[END_REF] for a general presentation and [START_REF] Van | Reduced order dynamical modelling of reaction systems: a singular perturbation approach[END_REF][START_REF] Hans | Analysis of Complex Reaction Networks in Signal Transduction[END_REF] for more chemically oriented texts. Our contribution [START_REF] Boulier | Model Reduction of Chemical Reaction Systems using Elimination[END_REF] consisted in formulating it, in the context of chemical reaction systems, as the result of a differential elimination process [START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF]. In this section, we present it over System [START_REF] Bächler | Thomas Decomposition of Algebraic and Differential Systems[END_REF], by showing how to obtain the famous Henri, Michaelis and Menten formula [START_REF] Victor | Lois générales de l'Action des Diastases[END_REF][START_REF] Michaelis | Die kinetik der invertinwirkung[END_REF]

dS dt (t) = - V max S(t) K + S(t) (4)
where V max and K are two parameters which can be expressed from the kinetic constants of the chemical reaction system. See [START_REF] Lena | Quasistationarität und Fast-Invariante Mengen Gewönlicher Differentialgleichungen[END_REF] for a closely related work.

Our algorithm was successfully applied over a more complicated model [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF][START_REF] Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF], featuring two genes, arising from the modeling problem of the circadian clock of the green algae Ostreococcus tauri [START_REF] Derelle | Genome Analysis of the smallest free-living eukaryote Ostreococcus tauri unveils many unique features[END_REF].

Computation using DifferentialAlgebra. Our method relies on the assumption that the chemical reactions are split into two sets: the set of the fast reactions, and the set of the slow ones. In the case of the Henri, Michaelis and Menten reduction, the revertible reaction is supposed to be fast, compared to the third one:

k 1 , k -1 ≫ k 2 .
One starts by building the natural deterministic model.

> with (LinearAlgebra): > X := <E(t), S(t), ES(t), P(t)>:

> V := <k[1]*E(t)*S(t), k[-1]*ES(t), k[2]*ES(t)>: > N := <<-1, -1, 1, 0> | <1, 1, -1, 0> | <1, 0, -1, 1>>: > X, N, V; [E(t) ] [-1 1 1] [ ] [ ] [k[1] E(t) S(t)] [S(t) ] [-1 1 0] [ ] [ ], [ ], [ k[-1] ES(t) ] [ES(t)] [ 1 -1 -1] [ ] [ ] [ ] [ k[2] ES(t) ] [P(t) ] [ 0 0 1]
Then, in the vector V , one replaces the laws corresponding to the fast reactions by two unknown laws F 1 (t) and F -1 (t). Denote W this new vector. The approximated system is obtained by enlarging the ODE system dX/dt = N • W , with the equation

k 1 E(t) S(t) = k -1 E(t) S(t)
, which gives the algebraic variety where the fast reactions would be equilibrated, if the slow reactions did not exist.

[--E(t) = -F[1](t) + F[-1](t) + k[2] ES(t), dt d d --S(t) = -F[1](t) + F[-1](t), --ES(t) = F[1](t) -F[-1](t) -k[2] ES(t), dt dt d --P(t) = k[2] ES(t), 0 = k[1] E(t) S(t) -k[-1] ES(t)] dt
This differential-algebraic equation can now be simplified, by means of a differential elimination software, such as the DifferentialAlgebra package or the recent [START_REF] Bächler | Thomas Decomposition of Algebraic and Differential Systems[END_REF]. With the DifferentialAlgebra package, the output involves three cases. We only give the equation which expresses the evolution of S(t), in the general case.

> with (DifferentialAlgebra): -------------------------------------------------

> Field := field (generators = [k[1],k[-1],k[2]]); Field := field(generators = [k[1], k[-1], k[2]]) > R := DifferentialRing (blocks = [[F[1],F[-1]], [ES,E,P,S], [k[1](),k[-1](),k[2]()]], derivations = [t]); R := differential_ring > ideal := RosenfeldGroebner (redsys, basefield = Field, R): ideal := [regular_differential_chain, regular_differential_chain, regular_differential_chain] > Equations (ideal[1], solved, leader=diff(S(t),t)); 2 2 d E(t) S(t) k[1] k[2] + E(t) S(t) k[1] k[-1] k[2] [--S(t) = -
] dt 2 E(t) k[1] k[-1] + S(t) k[1] k[-1] + k[-1]
This equation is not yet the Henri, Michaelis and Menten formula, since some minor hypotheses are missing: one still needs to take conservation laws (3) and some initial values into account. Last, one needs to rename K = k 1 /k -1 and V max = k 2 E(0) and to neglect some monomial K E(0), assuming S(0) ≫ E(0).

Here is a sequence of computations, which takes these hypotheses into account and leads to the sought formula. ---------] dt S(t) + K Computation using MABSys. The same computation can also be performed by the dedicated MABSys package [START_REF] Lemaire | The MABSys MAPLE package[END_REF][START_REF] Lemaire | Mabsys: Modeling and analysis of biological systems[END_REF], which relies on the MAPLE Regu-larChains [START_REF] Lemaire | The RegularChains library in MAPLE 10[END_REF] and ELPS [START_REF] Sedoglavic | Expanded Lie Point Symmetry (MAPLE package)[END_REF] packages. The MABSys package gathers as input a chemical reaction system. Each reaction can be defined as fast or slow. The package permits to obtain directly the stoichiometry matrix, the vector of the laws, the system of the equilibria, the conservation laws that can be read from the stoichiometry matrix and the natural deterministic model.

d sys := [--E(t) = -k1 E(t) S(t) + km1 C(t) + k2 C(t), dt d --S(t) = -k1 E(t) S(t) + km1 C(t), dt d d --C(t) = k1 E(t) S(t) -km1 C(t) -k2 C(t), --P(t) = k2 C(t)] dt dt
The package provides a function that performs the quasi-steady state approximation using the information provided by the fast boolean. -------------------------------------------- It also provides a function that permits to reduce the number of parameters of the system. Over this example, the parameter k 1 is removed, using the assumption that it is strictly positive. The change of coordinates is provided, following the syntax: "new parameter = rational function of the old parameters". ----------------------------------

- dt 2 2 S(t) + 2 km1 S(t) + E_0 km1 + km1 > output[2]; km1 [km1 = ---] k1
Other functions, which are still prototypes, hence not yet integrated in MABSys, permit to perform the final approximation which yields the Henri, Michaelis and Menten formula.

Reducing and Reparametrizing Models

This section, in which we give some more details on the method which permitted us to remove the k 1 parameter, is much inspired from [38, chapter 2] and [START_REF] Lemaire | A method for semi-rectifying algebraic and differential systems using scaling type Lie point symmetries with linear algebra[END_REF]. The goal consists in reducing and reparametrizing an initial model into a reduced one (the initial model possibly comes from the quasi-steady state approximation method). The reduced model is equivalent to the initial one, by an invertible change of coordinates, which preserves the most important properties: the positive steady points of the initial model are in bijection with the positive steady points of the reduced one, the initial model presents an oscillating behaviour if, and only if, the reduced model does, and so on. The computed changes of coordinates are given by Lie symmetries of the initial system of the simplest type: they are restricted to scalings. Moreover, the method distinguishes the parameters which are supposed to be strictly positive (dividing by them is allowed) from the ones which are only supposed to be nonnegative (dividing by them is forbidden). The computations of the scalings is performed by the ELPS package. It is important to point out that this package computes, in general, a restricted set of the scalings of the input system, in order to preserve a worst case polynomial complexity [START_REF] Sedoglavic | Reduction of Algebraic Parametric Systems by Rectification of their Affine Expanded Lie Symmetries[END_REF].

The method is illustrated over System [START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF], whose natural deterministic model [START_REF] Murray | Mathematical Biology I. An Introduction[END_REF]Eqs. (7.20)] depends on two variables and four parameters.

> ODS := [ > diff(A(t),t) = a -k1 * A(t) + k2 * A(t)^2 * B(t), > diff(B(t),t) = b -k2*A(t)^2*B(t) > ]; d 2 d 2 ODS := [--A(t) = a -k1 A(t) + k2 A(t) B(t), --B(t) = b -k2 A(t) B(t)] dt dt
The four parameters are supposed to be strictly positive. The steady points of the natural deterministic model are the zeros of non differential system obtained by equating to zero, the right hand sides of the model equations. The steady points obviously depend on the four parameters.

The reduction step. It makes use of the scalings of the natural deterministic model. These scalings can be described by the two following changes of coordinates. The first one depends on a parameter ν 1 . The second one depends on a parameter ν 2 .

(t, A, B, k 1 , k 2 , a, b) -→ t, A ν 1 , B ν 1 , k 1 , k 2 ν 2 1 , a ν 1 , b ν 1 . (t, A, B, k 1 , k 2 , a, b) -→ t ν 2 2 , A ν 2 , B ν 2 , k 1 ν 2 2 , k 2 , a ν 3 2 , b ν 3 2 .
Assuming that the four parameters are strictly positive, these scalings permit to rewrite the natural deterministic model in a new coordinate set. The new system depends on two parameters instead of four. The change of coordinates is given using the syntax: "new variables = rational fractions of the old ones".

> out := InvariantizeByScalings (ODS, [a,b,k1,k2]): > ODS_reduced := out [START_REF] Bächler | Thomas Decomposition of Algebraic and Differential Systems[END_REF];

ODS_reduced := d 2 d 2 [--A(t) = 1 -A(t) + k2 A(t) B(t), --B(t) = b -k2 A(t) B(t)] dt dt > Change_of_coordinates := out[2]; 2 k2 a A k1 B k1 Change_of_coordinates := [b = b/a, k2 = -----, t = t k1, A = ----, B = ----] 3 a a k1
The reparametrization step. It makes use of the scalings of the system which gives the steady points of ODS reduced. These scalings can be described by the following change of coordinates, which depends on a parameter ν 3 .

(A, B, k 2 , b) -→ A, B ν 3 , k 2 ν 3 , b .
The key idea is now very simple: apply these scalings, not on the system which gives the steady points of ODS reduced, but on ODS reduced itself. In the resulting system, the parameter k 2 is still present. However, it has no control anymore on the location of the steady points: it can only control their stabilities.

> out := SemiRectifySteadyPoints (ODS_reduced, [b,k2,A,B], []): > Reparamatrized_system := out[1][1]; Reparamatrized_system := d 2 d 2 [--A(t) = 1 -A(t) + A(t) B(t), --B(t) = (b -A(t) B(t)) k2] dt dt > Change_of_coordinates := out[1][3]; Change_of_coordinates := [B = k2 B]
The dynamics of chemical reaction systems can also be studied with a stochastic point of view. The parameters of the reactions are stochastic constants, which are probabilities that reactions occur per unit of time [17, page 2342]. To each species A, one associates a random variable A(t) which counts the number of molecules of A. In a celebrated paper [START_REF] Daniel | Exact Stochastic Simulation of Coupled Chemical Reactions[END_REF], Gillespie provided a numerical simulation algorithm, based on a strong rigorous analysis, under some assumptions ; the main one being that the chemically reacting system keeps being well-stirred over the time.

In the context of biological modeling, Gillespie's main assumption is very unlikely. However, the analysis of the stochastic behaviour of chemical reaction systems is a very important counterpart of the deterministic analysis. In a living cell, some molecules may occur in a very little number: if one views a gene, in the "active" state, as a chemical species, the number of molecules is zero or one, and handling this number as a real valued concentration is contestable. In such cases, it is known that there exist many different stochastic behaviours which are, at least, difficult to reproduce by a deterministic model. See [19, page 454]. Moreover, stochastic simulations may reproduce the surprising effects of the "noise" on the models, as pointed out in [START_REF] José | Mechanisms of noise-resistance in genetic oscillators[END_REF] and the references therein.

It is well-known that the evolution, over the time t, of the moments (mean, variance, covariance) of the random variable associated to each species may be described by a system of ordinary differential equations, at least for chemical systems the reactions of which do not involve more than one reactant [START_REF] Paulsson | Models of stochastic gene expression[END_REF][START_REF] Paulsson | Stochastic Modeling of Intracellular Kinetics[END_REF]. For more general systems, the differential system is, in general, infinite, and approximating it by a finite ODE system, by performing a so-called moment closure, is a difficult problem [START_REF] Singh | Lognormal Moment Closures for Biochemical Reactions[END_REF][START_REF] Gillespie | Moment-closure approximations for mass-action models[END_REF]. In this area, our contribution [START_REF] Vidal | Models of Stochastic Gene Expression and Weyl Algebra[END_REF] consisted in showing how to build these ODE systems by using Weyl algebra methods. The use of Weyl algebra led us to a few algorithmic improvements (reducing expressions swells and taking advantage of linear conservation laws at an early stage of the formula generation process).

Analysis of Statistical Moments

Before proceeding, one needs to associate a precise stochastic dynamics to a chemical reaction system. This can be achieved using stochastic Petri nets1 , endowed with their standard temporisation [START_REF] Vidal | Models of Stochastic Gene Expression and Weyl Algebra[END_REF] ; or by relying on the Gillespie's stochastic simulation algorithm.

Define the state ν of a given chemical reaction system as a vector of nonnegative integers (one coordinate per chemical species, each coordinate being a number of molecules) ; π ν (t) as the probability that the chemical reaction system be in the state ν at time t ; and the probability generating function

φ(t, z) = ν≥0 π ν (t) z ν , (5)
where z ν stands for z ν1 1 z ν2 2 • • • z νn n (one "symbol" z i per chemical species). Given any chemical reaction system, it is possible to compute a general equation for φ [33, Eq. (5.60)] as follows. To each reaction

α 1 A 1 + α 2 A 2 + • • • + α n A n c ----→ β 1 A 1 + β 2 A 2 + • • • + β n A n , associate the differential operator c α! z β -z α ∂ ∂ z α • (6)
The very same operator can be denoted using Euler operators

θ i = z i ∂/∂ z i instead of partial derivatives c α! z β-α -1 θ α , (7)
where θ α denotes the product of the θ

α i i = θ i (θ i -1) • • • (θ i -α i +1)
. Define H as the sum of the differential operators [START_REF] Boulier | On proving the absence of oscillations in models of genetic circuits[END_REF], or [START_REF] Brown | QEPCAD B: a program for computing with semialgebraic sets using CADs[END_REF], for all reactions of the considered chemical reaction system. The general equation for φ is:

∂ ∂ t φ(t, z) = H φ(t, z) . (8)
With partial derivatives, the differential operator H, for System (2), is

a (z A -1) + k 1 (1 -z A ) ∂ ∂ z A + k 2 2 z 3 A -z 2 A z B ∂ 3 ∂ z 2 A ∂ z B + b (z B -1) .
With Euler operators, the differential operator H, for System (1), is

k 1 z ES z E z S -1 θ E θ S + k -1 z E z S z ES -1 θ ES + k 2 z E z P z ES -1 θ ES .
By differentiating the probability generating function [START_REF] Boulier | Applying a rigorous quasi-steady state approximation method for proving the absence of oscillations in models of genetic circuits[END_REF] and evaluating it at

z 1 = z 2 = • • • = z n = 1
, one gets formulae which bind φ and the moments of the random variables which count the numbers of molecules. A mere evaluation yields: φ(t, z) | z=1 = 1. Differentiating (5) with respect to any z i and evaluating at z = 1 provides the expected value of the number of molecules of the chemical species associated to z i . These ideas are illustrated using a prototype software available at [START_REF] Petitot | The MAGNUS MAPLE software[END_REF], on System [START_REF] Bächler | Thomas Decomposition of Algebraic and Differential Systems[END_REF]. Variables are numbered, with the convention (E, S, ES , P ) = (1, 4, 2, 3). The differential operator H is displayed using Euler operators.

/ z[2] \ /z[1] z[4] \ H := k[1] |----------1| theta[1] theta[4] + k[-1] |----------1| theta[2] \z[1] z[4] / \ z[2] / /z[1] z[3] \ + k[2] |----------1| theta[2] \ z[2] /
The next command permits to generate an ODE system for the expected values and the covariances of the random variables which count the molecules of E, S, ES and P (the variable x 1 (t) is the expected value of the random variable which counts the molecules of the enzyme E ; the variable x 1,4 (t) is the covariance of the random variables which count, respectively, the molecules of the enzyme E and the product P ). Since some reactions of System (1) involve two reactants, the ODE system is infinite. The function truncates it. The parameter iv receives some information on initial values, for a further numerical integration.

> ODEs := mean_equations (HSyst, 'iv');

Conclusion

Beyond the improvements of the underlying theories, a lot of improvements can be brought to the software. Since this paper is software oriented, let us focus on that single issue. First, our packages could certainly be polished and merged. Second, they miss very important tools dedicated to the solving problem of polynomial systems in the field of the real numbers, such as [START_REF] Dolzmann | Redlog: computer algebra meets computer logic[END_REF][START_REF] Brown | QEPCAD B: a program for computing with semialgebraic sets using CADs[END_REF][START_REF] Xia | DISCOVERER: A tool for solving semi-algebraic systems[END_REF][START_REF] Safey | RAGLib: A library for real solving polynomial systems of equations and inequalities[END_REF][START_REF] Moroz | DV: Un logiciel de classification des solutions réelles d'un système paramétré[END_REF]. See [START_REF] Niu | Qualitative Analysis of Biological Systems Using Algebraic Methods[END_REF] for a recent study of the sizes of the problems that can be investigated using some of the best available such tools. However, even the nice integrated MAPLE package that we could foresee would miss two features: sort of a model checking functionality, which would permit to the practitioner to query models, as [START_REF] Fages | Modelling and querying interaction networks in the biochemical abstract machine BIOCHAM[END_REF][START_REF] Batt | Efficient parameter search for qualitative models of regulatory networks using symbolic model checking[END_REF] do ; and a user interface of the same quality as that of, say, Cytoscape [START_REF]Cytoscape: An Open Source Platform for Complex Network Analysis and Visualization[END_REF]. These two missing features, which are both user interface related, may seem a bit irrelevant to the traditional audience of a computer algebra conference.

The need for them is however pretty obvious for any lecturer, teaching symbolic techniques close to ours, among other systems biology approaches.

  > W := <F[1](t), F[-1](t), k[2]*ES(t)>: > equilibre := 0 = k[1]*E(t)*S(t) -k[-1]*ES(t): > redsys := map (diff, X, t) = N . W: > redsys := [ seq (lhs (redsys) [i] = rhs (redsys) [i], i = 1 .. Dimension (X)), equilibre ]; d redsys :=

>

  conservation_laws := [E(t) + ES(t) = E0 + ES0, S(t) + ES(t) + P(t) = S0 + ES0 + P0]: > hypotheses := [P0 = 0, ES0 = 0, op (conservation_laws)]: > R := DifferentialRing ( blocks = [[F[1],F[-1]], [ES,E,P,S], [ES0(),E0(),P0(),S0()], [k[1](),k[-1](),k[2]()]], derivations = [t]): > ideal := RosenfeldGroebner ([ op(redsys), op(hypotheses) ], R, basefield = Field): > formula := Equations (ideal[1], solved, leader = diff (S(t),t)): > formula := subs (k[-1]=K*k[1], formula): > formula := normal (formula): > formula := algsubs (k[2]*E0=Vmax, formula): > formula := normal (subs (K*E0=0, formula)); d Vmax S(t) formula := [--S(t) = -

>

  with(MABSys): > R1 := NewReaction (E+S, C, MassActionLaw(k1), fast=true): > R2 := NewReaction (C, E+S, MassActionLaw(km1), fast=true): > R3 := NewReaction (C, E+P, MassActionLaw(k2)): > RS := [R1,R2,R3]: > ConservationLaws(RS); [C + P + S -C_0 -P_0 -S_0, C + E -C_0 -E_0] > sys := ReactionSystem2ODEs(RS, [E,S,C,P]);

>

  output := ModelReduce(RS, [E,C,P,S], useConservationLaws=true): > red_sys := output[1][1]: > red_sys := subs (C_0=0, P_0=0, red_sys): > red_sys[4]; d E_0 k2 k1 S(t) (k1 S(t) + km1) --S(t) = -

>

  output := InvariantizeByScalings (red_sys, [k1], [km1,k2,E,C,P,S]): > red_sys2 := output[1]: > red_sys2[4]; d S(t) k2 E_0 (S(t) + km1) --S(t) = -

>

  cr1 := ChemicalReaction (E+S = ES, k[1], ""): > cr2 := ChemicalReaction (ES = E+S, k[-1], ""): > cr3 := ChemicalReaction (ES = E+P, k[2], ""): > CRS := ChemicalSystem (cr1,cr2,cr3): > HSyst := HamiltonianSystem (CRS, z, theta, x(t)); HSyst := hamiltonian_system > map (numero, [E,S,ES,P], HSyst); [1, 4, 2, 3] > H := Hamiltonian (HSyst);

d

  ODEs := [--x[1](t) = -%2 + (k[-1] + k[2]) x[2](t) -%1, dt d d --x[2](t) = %2 + (-k[-1] -k[2]) x[2](t) + %1, --x[3](t) = k[2] x[2](t), dt dt d d --x[4](t) = -%2 + k[-1] x[2](t) -%1, --x[1, 1](t) = ... dt dt %1 := k[1] x[1, 4](t) %2 := k[1] x[4](t) x[1](t)

According to some sources, cited on the Wikipedia, Petri nets were invented by Carl Adam Petri to model chemical reaction systems at the age of 13.

⋆ This work has benefited from the support of the French ANR (decision number ANR-2010-BLAN-0109-03).


We consider a compressible fluid occupying an open bounded set depending on time.

∂ρu ∂t

+ div (ρu ⊗ u) -µ∆u + a∇ρ γ -ξ∇div u = ρf , on Q (1)

∂ρ ∂t + div (ρu) = 0 on Q (2)
where µ > 0, µ + ξ > 0, a > 0, T > 0, γ ∈ (1, ∞).

Boundary condition :

u(X + d(X , t), t) = ∂d ∂t on Σ
where d is a given smooth function describing the Σ evolution.

Some preliminary remarks :

On boundary conditions

Usually, for non homogeneous problem, it is necessary to impose value of density when u • n < 0. But, in our problem, even if we have non homogeneous boundary condition on the velocity, no fluid go across the boundary and it is no necessary to impose the value of density when u • n < 0.

On Eulerian or Lagrangian description

The more natural method to treat the difficulty linked to a moving domain is to use a lagrangian description of the equations. But we are not able to obtain energy estimate using this Lagrangian form.

The homogeneous problem

As usual for a non homogeneous Dirichlet problem, we assume that there exists a sufficiently smooth function w on Q such that w = ∂d ∂t on γt and verifying

w ∈ L 2 (0, T , H 1 (Ωt )) ∩ L ∞ (0, T ; L 2γ/(γ-1) (Ωt )) div w ∈ L 2 (0, T ; L ∞ (Ωt ))
∂w ∂t in a good space.

By tacking v = uw , we can prove the following a priori result :

Lemma

Let (v , ρ) a classical solution of the problem, then under some assumption on initial data and for a finite time T , v and ρ verify :

v ∈ L 2 (0, T ; H 1 0 (Ωt )), ρ ∈ L ∞ (0, T ; L γ (Ωt )), √ ρv ∈ L ∞ (0, T ; L 2 (Ωt )).

Proof

A large part of the estimate is usual and we give only a sketch of proof to treat the pressure term av ∇ρ γ . To handle this one, we multiply the mass equation by a γ γ-1 ρ γ-1 :

a γ -1 ∂ρ γ ∂t + γa γ -1 ρ γ-1 div (ρu) = 0, with div (ρρ γ-1 u) = ρ γ-1 div (ρu) + γ -1 γ u∇ρ γ .
After some straight-forward manipulations and an integration on Ωt , we obtain

Z Ωt av ∇ρ γ dx = a Z Ωt ρ γ div w dx + a γ -1 Z Ωt ∂ρ γ ∂t dx + a γ -1 Z γt ρ γ w • n dS.
Finally, by using the Leibniz formula, the boundary terms vanish, and we obtain the following result :

Z Ωt av ∇ρ γ dx = a Z Ωt ρ γ div w dx + a γ -1 d dt Z Ωt ρ γ dx.
Energy estimate is the following

d dt Z Ωt ρ |v | 2 2 + a γ -1 ρ γ dx + Z Ωt µ|Dv 2 |dx = Z Ωt ρ " f - ∂w ∂t « v + (ρu • ∇)w • v -aρ γ div w + µdiv (v ∇w ) -µDvDw dx.
The first two terms of the right-hand side are easily estimated by the Cauchy Schwarz inequality For the two remaining terms, we use the bound on div w in L 1 (0, T ; L ∞ (Ωt )) and we can write

Z Ωt ρ γ div wdx ||div w || L ∞ (Ωt ) Z Ωt ρ γ dx, Z Ωt (ρu • ∇)w • vdx C ||div w || L ∞ (Ωt ) »Z Ωt ρv 2 dx + ||ρ|| L γ (Ωt ) ||w || 1 2 L 2γ/(γ-1) (Ωt )
-.

Then we the deduce from Grönwall's lemma the desired bounds.

Time discretization

Proof of the existence is obtained by building a sequence of approached solutions based on the time discretisation of the continuous equations. The main difficulty is linked to the time depending domain. A classical Eulerian discretization requires to evaluate the value at a given position x i at time t k depending of the value at the same position at time t k-1 . But this point can be out of the domain Q.

∂ρ ∂t (x i , t k ) ≃ ρ(x i , t k ) -ρ(x i , t k-1 ) ∆t
For M 1, we denote ∆t = T M and we define for k 0 a sequence of (m + 1) domains Ω k by taking

Ω k = Q ∩ k∆t. Let ρ k defined on Ω k , we define ρk-1 a function also defined on Ω k such that ρk-1 = ρ k-1 on (Ω k-1 ∩ Ω k ) and ρk-1 = 0 on (∁ Ω k (Ω k-1 ∩Ω k ) ).
Similarly, we define ūk-1 a function such that ūk-

1 = u k-1 on (Ω k-1 ∩ Ω k ), but ūk-1 are non imposed on (∁ Ω k (Ω k-1 ∩Ω k ) )
. Since this term is systematically multiplied par ρk-1 , the value of the product is zero in (∁

Ω k (Ω k-1 ∩Ω k )
) for all chosen value of ūk-1 Let w k = w (k∆t) where w is the function defined previously, by induction for k 1, we define (ρ k , v k ) to be the solution of our discretized problem on Ω k :

ρ k -ρk-1 + ∆tdiv (ρ k u k ) -C k-1 = 0 ρ k v k -ρk-1 vk-1 + ρk-1 (w k -wk-1 ) + w k (ρ k -ρk-1 ) + ∆tdiv (ρ k u k ⊗ u k ) -µ∆t∆u k + ∆t a∇ρ γ k = 0 where C k-1 = 1 meas Ω k Z ∁ Ω k-1 (Ω k-1 ∩Ω k ) ρ k-1 dx and u k = v k + w k .
At this time, we have to considere existence of solution for this stationary problem. This one can be obtained by a fixed point argument after obtention of some compactness results.

Compactness results

Lemma Let (ρ k , v k ) (0 k M) a sequence of weak solutions of our problem with the homogeneous boundary conditions on v k and some initial conditions ρ 0 0, v 0 a.e. in Ω 0 . Let

ρ k = 0 on {x ∈ ∂Ω k , w (x) • n < 0} sup 0 l M ||ρ l || L γ (Ω l ) + ||ρ l |v l | 2 || L 1 (Ω l ) C M X k=1 ∆t||∇v l || 2 L 2 (Ω l ) C M X k=1 || p ρk-1 (v k -vk-1 )|| 2 L 2 (Ω k ) C M X k=1 Z Ω k [(γ -1)ρ γ k + ργ k-1 -γρ γ-1 k ρk-1 ]dx C M X k=1 Z Ω k |ρ k -ρk-1 | γ dx C
Where, here and below, C denotes various positive constants independent of k, ∆t and thus M.

Proof

We give only one step of the demonstration to trait the pressure term.

To do this, we multiply the discretized mass equation by a ρ γ-1 k γ -1 and after some manipulations, we deduce :

a∆tv k ∇ρ γ k = γa γ -1 ρ γ k - γa γ -1 ρ γ-1 k ρk-1 + a∆T γ γ -1 div (u k ρ γ k ) -a∆tw k ∇ρ γ k - γC k-1 γ -1 ρ γ-1 k .
Then, after integration over Ω k , we obtain a∆t

Z Ω k v k ∇ρ γ k dx = a γ -1 Z Ω k ρ γ k dx - a γ -1 Z Ω k-1 ρ γ k-1 dx + a γ -1 Z ∁ Ω k-1 (Ω k-1 ∩Ω k ) ρ γ k-1 dx + aγ∆T γ -1 Z Ω k div (u k ρ γ k )dx -a∆T Z Ω k w k ∇ρ γ k dx - γC k-1 γ -1 Z Ω k ρ γ-1 k dx + a γ -1 Z Ω k (γ -1)ρ γ k + ργ k-1 -γρ γ-1 k ρk-1 dx By adding from k = 1 to M we obtain 1 2 || √ ρ M v M || 2 L 2 (Ω M ) + 1 2 M X k=1 Z ∁ Ω k-1 (Ω k-1 ∩Ω k ) ρ k-1 v 2 k-1 + 3µ 2 M X k=1 ∆t||∇v k || 2 L 2 (Ω k ) + a γ -1 ||ρ M || γ L γ (Ω M ) + " a γ -1 -ε « M X k=1 ||ρ k-1 | γ L γ (∁ Ω k-1 Ω k-1 ∩Ω k ) + M X k=1 Z Ω k [(γ -1)ρ γ k + ργ k-1 -γρ γ-1 k ρk-1 ] + 1 2 M X k=1 || (1 -∆t) p ρk-1 (v k -vk-1 )|| 2 L 2 (Ω k ) 1 2 || √ ρ 0 v 0 || 2 L 2 (Ω 0 ) + a γ -1 ||ρ 0 || γ L γ (Ω 0 ) + M X k=1 ∆t » 3 2 ||Dw k || L ∞ (Ω k ) + ||f k || 2 L 2γ/(γ-1) (Ω k ) || - || √ ρ k v k || 2 L 2 (Ω k ) + M X k=1 ∆t » ||Dw k || L ∞ (Ω k ) 1 2 ||w k || 1 2 L 2γ/(γ-1) (Ω k ) -aC ||Dw k || L ∞ (Ω k ) + 3 2 ||f k || 2 L 2γ/(γ-1) (Ω k ) + Cε - ||ρ k || L γ (Ω k ) -µ M X k=1 ∆t 1 2 ||∇w k || 2 L 2 (Ω k ) + 3 2 ||f 0 || 2 L 2γ/(γ-1) (Ω 0 ) ||ρ 0 || L γ (Ω 0 )
Using the discrete Gronwall lemma, we deduce the bound on ρ k v 2 k , ρ k and ∇v k .

To ensure convergence of the solution of our equations to this of the continuous problem, we need to obtain an other bound for ρ to treat the pressure term. We now can prove the following lemma

Lemma If γ > 2 for N = 2 or γ ≥ 3 for N = 3, we have ∀ϕ ∈ C ∞ ([0, T ]; C ∞ 0 (Ωt )) M X k=1 ∆t||ϕρ k || γ+1 L γ+1 ≤ C

Proof

The proof follow these proposed by P.L. Lions by applying the ∆ -1 div operator at the discrete momentum equation. To apply this operator in case of Dirichlet boundary condition, we need to localize the equations by applying a cut-off function. In our case of moving domain, we define for any fixed compact set

K ⊂ Q, a cut-off function ϕ ∈ C ∞ ([0, T ]; C ∞ 0 (Ωt )
) such that ϕ ≡ 1 on K and ϕ 0 on Q. The end of the proof corresponds to this given by P.L. Lions.

Passage to the limit

Introducting for any function (ϕ 0 , ϕ 1 , ...ϕ M ) the piecewise constant (respectively piecewise affine) function φ and φ (respectively φ) defined by

φ(x, t) = ϕ k+1 (x), if k∆t t < (k + 1)∆t φ(x, t) = φk (x), if k∆T t < (k + 1)∆t (respectively φ(x, t) = ϕ k (x) + (t -k∆t) " ϕ k+1 -φk ∆t « (x), if k∆t t < (k + 1)∆t)
We have

Z T 0 φ(t) dt = M X k=1 ∆t ϕ k , Z T 0 φ(t) dt = M X k=1 ∆t ϕ k-1 + ϕ k 2 and ∂ φ ∂t = ϕ k+1 -φk ∆t if k∆t t < (k + 1)∆t
With this notation, we can observe that our discretized equations on [(k -1)∆t, k∆t[, (k 1)

ρ k -ρk-1 + ∆tdiv (ρ k u k ) -C k-1 = 0 ρ k v k -ρk-1 vk-1 + ρk-1 (w k -wk-1 ) + w k (ρ k -ρk-1 ) + ∆tdiv (ρ k u k ⊗ u k ) -µ∆t∆u k + ∆t a∇ρ γ k = 0
can by written as

∂ ρ ∂t + div (ρû) = C ∆t ∂f ρu ∂t + div (ρû ⊗ û) -µ∆û + a∇ρ γ = 0 where û = v + ŵ and f ρu = f ρv + f ρw .
It is easy to see that the term C ∆t goes to zero when ∆t goes to zero. We recall that

C k-1 = 1 meas Ω k Z ∁ Ω k-1 (Ω k-1 ∩Ω k ) ρ k-1 dx. Z ∁ Ω k-1 (Ω k ∩Ω k +1) ρ k-1 dx ||ρ k-1 || L γ (∁ Ω k-1 (Ω k-1 ∩Ω k ) ) meas (∁ Ω k-1 (Ω k ∩Ω k +1) ) γ γ-1
We introduce then the d k sequence such that d 0 = 0 and for all k and all point X belonging on the boundary of

Ω k , d k (X ) = d k-1 (X ) + R k∆t (k-1)∆t w (X , t)dt. Then meas (∁ Ω k-1 (Ω k ∩Ω k -1) ) C Z γ 0 |d k -d k-1 |dS C ′ ∆T
where C and C ′ depend of the regularity of the domain.

Then C ∆t (C ′ ∆T ) 1 γ-1 ||ρ k-1 || L γ (∁ Ω k-1 (Ω k-1 ∩Ω k ) ) C ′′ ∆T 1 γ-1 From the estimate sup 0 l M ||ρ l || L γ (Ω l ) C and M X k=0 ∆t||ρ k || γ+1
L γ+1 ≤ C , we may assume that, extracting subsequences if necessary, as ∆t goes to 0, ρ converge weakly in L ∞ (0, T ; L γ (Ωt )) (weak-*) and in L γ+1 (L γ+1 loc (Ωt )) for some ρ 0.

In addition, from

M X k=1 Z Ω k |ρ k -ρk-1 | γ dx C we deduce that Z T 0 Z Ωt |ρ(t) -ρ(t)| γ dxdt = ∆t M X k=1 Z Ω k |ρ k -ρk-1 | γ dx C ∆t
and we see easily that ρ -ρ converge to 0 as ∆t goes to 0+ in L q (0, T ; L γ ) for all 1 q < ∞. Moreover,

ρ -ρ = ρ k -ρk-1 ∆t (t -(k -1)∆t) | {z } ∆t on [(k -1)∆t; k∆t[, then, Z T 0 Z Ωt |ρ -ρ| γ dxdt Z T 0 Z Ωt |ρ k -ρk-1 | γ dxdt ||ρ -ρ|| γ L γ (Q)
and we deduce that ρ -ρ converge to 0 as ∆t goes to 0+ in L q (0, T ; L γ ) for all 1 q < ∞.

Passage to the limit in non linear terms

To ensure convergence of non linear terms, we recall the following general lemma (P.L. Lions)

Lemma

Let gn, hn converge weakly to g , h respectively in L p1 (O, T ; L p2 (Ω)), L q1 (O, T ; L q2 (Ω)) where 1 p1, p2 ∞,

1 p1 + 1 q1 = 1 p2 + 1 q2 = 1.
We assume in addition that ∂gn ∂t is bounded in L 1 (0, T ; W -m,1 (Ω)) for some m 0 independent of n. ||hnhn(. + ξ, t)|| L q1 (O,T ;L q2 (Ω)) → 0 as |ξ| → 0, uniformly in n.

Then, gnhn converge to gh (in the sens of distribution in Ω × (0, T )).

The estimate sup

0 l M ||ρ l || L γ (Ω l ) + ||ρ l |u l | 2 || L 1 (Ω l ) C imply that √ ρû, g √ ρu are bounded in L ∞ (0, T ; L 1 (Ωt )) and √ ρ, √ ρ are bounded in L ∞ (0, T ; L 2γ (Ωt )). ⇒ ρû, f ρu are bounded in L ∞ (0, T ; L 2γ/(γ+1) (Ωt )).
We can also claim that ρû -ρū and ρûf ρu converge to 0 as ∆t goes to 0+ in L q (0, T ; L r (Ωt )) with 2 q < ∞, 1 r < 2γ γ+1 . This convergence is obvious using previous bounds and observing that ρ

k u k -ρk-1 ūk-1 = (ρ k -ρk )u k + ρk-1 (u k -u k-1 ) From ∂ ρ ∂t + div (ρû) = C
∆t and using the bound on ρû, we deduce that ∂ ρ ∂t is bounded in L ∞ (0, T ; W -1,2γ/(γ-1) (Ωt )) ⊂ L 1 (0, T ; W -1,1 ). Then, applying previous lemma we deduce that ρû converge to ρu (in the sense of distributions on Q).

This equation can then be written as follow :

∂ ρ ∂t + div (ρû) -div ((ρ -ρ)û) = C
∆t and, since (ρ -ρ) goes to 0 we can pass to the limit in this equation in the sense of distribution on Q to obtain ∂ρ ∂t + div (ρu) = 0

To pass to the limit in momentum equation, we first observe that |ρû ⊗ ûf ρu ⊗ û| |û||ρû -ρū|

We deduce then that ρû ⊗ ûf ρu ⊗ û converge to 0 as ∆t goes to 0+ in L 1 (0, T ; L r (Ωt )) ∩ L q (0, T ; L 1 (Ωt )) for all 1 r r 1 if N 3, 1 r < γ if N = 2 and for all q ∈ [1, +∞ with 1

r 1 = 1 γ + N-2 N
In view of bounds show above, we may prove that ∂ f ρu ∂t is bounded in L ∞ (0, T ; W -1,1 (Ωt )) + L 2 (0, T ; H -1 (Ωt )) ⊂ L 1 (0, T ; W -1,1 (Ωt )). Using the Lemma 5.1 (PL2), we deduce that f ρu ⊗ û, and thus ρ(û ⊗ û) converge weakly in L 1 (0, T ; L r (Ωt )) ∩ L q (0, T ; L 1 (Ωt )) to ρu ⊗ u.

To finish the proof, we need then to obtain a strong convergence of ρ in L γ (Q). We follow for thus the proof proposed by P.L. Lions in the case were ρ is bounded in L γ+1 (Q) ∩ L ∞ (0, T , L s (Ωt )) for some s > N, which is precisely the case here.

Then, when ∆t goes to 0+, our discretized system converge to the continuous system.

References :

For a fixed domain, mains result have be presented by P.L. Lions (Mathematical topics in fluid mechanics, volume 2, Compressible models, Oxford science Publication, 1996) For time depending domain, some works have been realized by B. Desjardins and M.J. Esteban, and by M. Boulakia when a fluid is contained in a fixed bounded set and some elastic structures are immersed in this fluid.

The problem of interaction between a compressible fluid bounded by a elastic plate is considered F. Flori and P. Orenga with the quasi Navier Stokes equations.

An existence result for a free boundary shallow water model is proposed by F. Flori, P. Orenga and M. Peybernes using a Lagrangian scheme.

Perspectives 1) This work is restricted at the case of given moving domain (Q thus d are given). We propose to extend this result at a problem where d is unknown, for instance, d can be a solution of a PDE describing movement of a plate or a shell.

The main difficulty is then that we need to construct the domain at each time step and the solution can be obtained by a fixed point method. Moreover, the equation verified by d need to give sufficient regularity to obtain a sufficient regularity on boundary Σ that permit to apply all previous result.

2) Developments of a numerical code based on the constructive method developed here.

1 How to calculate W ?

2 Since we do not pass numerically to the limit, how many manage ρ on the boundary.