Andrzej Mizera
email: amizera@ippt.gov.pl

Barbara Gambin
email: bgambin@ippt.gov.pl

Stochastic modelling of the eukaryotic heat shock response

Keywords: Mathematical modelling, Stochastic model, Computer simulations, Markov chain, Gillespie algorithm, Stationary distribution, Steady-state, Clustering

The heat shock response (HSR) is a highly evolutionarily conserved defence mechanism allowing the cell to promptly react to elevated temperature conditions and other forms of stress. It has been subject to intense research for at least two main reasons. First, it is considered a promising candidate for deciphering the engineering principles underlying regulatory networks.

Second, heat shock proteins (main actors of the HSR) play crucial role in many fundamental cellular processes. Therefore, profound understanding of the heat shock response would have far-reaching ramifications for the cell biology.

Recently, a new deterministic model of the eukaryotic heat shock response has been proposed in the literature. It is very attractive since it consists of only the minimum number of components required by any functional regulatory network, while yet being capable of biological validation. However,

Introduction

The heat shock response is the most highly evolutionarily conserved defence mechanism [START_REF] Lindquist | The heat-shock proteins[END_REF]. It exists in all eukaryotic cells, protects them from the damaging influence of elevated temperature and allows them to promptly react to other forms of environmental stress. The

A c c e p t e d m a n u s c r i p t

heat shock response has been subject to intense research recently, see [START_REF] Chen | Heat shock paradox and a new role of heat shock proteins and their receptors as antiinflammation targets[END_REF]Powers and Workman, 2007;[START_REF] Tijms | Chaperone regulation of the heat shock protein response[END_REF], for at least two reasons. On one hand, as a well-conserved mechanism, it is considered a promising candidate for deciphering the engineering principles underlying regulatory networks. On the other hand, heat shock proteins play crucial roles in many fundamental cellular processes such as protein biogenesis, dismantling of damaged proteins, activation of immune responses and signalling, see [START_REF] Kampinga | Thermotolerance in mammalian cells: protein denaturation and aggregation, and stress proteins[END_REF]Pockley, 2003). Therefore, profound understanding of the heat shock response would have far-reaching ramifications for the cell biology and could potentially allow for treatment of a number of diseases, such as neurodegenerative and cardiovascular disorders, cancer, ageing, see [START_REF] Balch | Adapting proteostasis for disease intervention[END_REF][START_REF] Liu | Overcomming immune toerance to A c c e p t e d m a n u s c r i p t cancer by heat shock protein vaccines[END_REF][START_REF] Lukacs | Heat shock proteins in cancer therapy[END_REF][START_REF] Morimoto | Proteotoxic stress and inducible chaperone networks in neurodegenerative disease and aging[END_REF][START_REF] Workman | Putting the heat on cancer[END_REF].

Although a number of mathematical models describing the heat shock response both in eukaryotes as well as in bacteria have been presented in the literature, see [START_REF] Donati | Oxidative injury and the heat shock response[END_REF]; [START_REF] Jones | Fast events in protein folding initiated by nanosecond laser photolysis[END_REF]; [START_REF] Parsell | The function of heat-shock proteins in stress tolerance: degradation and reactivation of damaged proteins[END_REF]; [START_REF] Peper | A mathematical model of the hsp70 regulation in the cell[END_REF]; Petre et al. (2009b); Remondini et al. (2006); Rieger et al. (2005); [START_REF] Szymańska | Mathematical modeling of heat shock protein synthesis in response to temperature change[END_REF], still a comprehensive mechanistic understanding of this process is lacking. In Petre et al. (2009b) a new model of the eukaryotic heat shock response together with an associated continuous mathematical model based on ordinary differential equations have been discussed. The novelty of the model in Petre et al. (2009b) is due to the fact that, unlike other previous models, it is based solely on welldocumented reactions and does not incorporate modelling "blackboxes" such as hypothetical, experimentally unsupported cellular mechanisms whose only

A c c e p t e d m a n u s c r i p t

purpose is to enforce appropriate behaviour. The simplified version of the model (see Petre et al. (2009a) for details) includes the temperature-induced protein misfolding, all three forms of heat shock factors: monomers, dimers and trimers, the backregulation of the transactivation of the heat shock protein encoding gene and the chaperone activity of heat shock proteins. At the same time, it contains as few reactions and reactants as possible. It is worth noticing that the model consists of only the minimum number of components required by any functional regulatory network: an activation mechanism and a feedback mechanism. Nonetheless, the associated continuous model predictions correlate well with experimental observations on the heat-induced transactivation of the hsp-encoding genes at different temperatures from the range 37 • C -43 • C (in particular, the prolonged transcription at 43 • C is confirmed) and the return to the original level of hsp production once the stress is removed (publication in preparation). Moreover, the model perfectly illustrates the experimentally observed process of "self-learning" of the HSR system: the response to a second consecutive heat shock is significantly lower.

This is due to a transient increase in the free hsp level caused by the preliminary heat shock. In other words, the increase is a form of temporary memory of the fact that the cell was recently exposed to heat shock conditions.

However, the undertaken modelling approach that utilises ordinary differential equations is just one of many other modelling paradigms (e.g. stochastic formulation, process calculi, Petri nets, etc.), which could be exploited in the context of the heat shock response. In this paper we follow one of the in Petre et al. (2009a). According to current scientific knowledge, ignoring quantum mechanical effects, biological systems can be viewed as deterministic of their very nature, with their dynamics entirely specified, given sufficient information on the state of the system (position, orientation and momentum of every single molecule) and a complete understanding of the chemistry and physics of the interactions between biomolecules. Unfortunately, we are still unable to model biological systems of realistic complexity and size using such a molecular dynamic approach [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF]). Therefore the current models admit far-reaching simplifications, which result in a higher level view of the system being modelled. However, these abstractions change the character of the dynamics, which becomes intrinsically stochastic and requires consideration of statistical physics to describe the stochastic process governing it. Especially at low concentrations of the involved reactants, random fluctuations may have a significant impact on the reaction dynamics, but the deterministic approach to chemical kinetics fails to capture such phenomena, see [START_REF] Mcadams | It's a noisy business! Genetic regulation at the nanomolar scale[END_REF]; [START_REF] Srivastava | Stochastic vs. deterministic modeling of intracellular viral kinetics[END_REF]. For example, let us consider the famous Lotka-Volterra system of coupled ordinary differential equations describing an ecological predator-prey model. The solutions of this system are known to be periodic (except for the stationary point) independently of the initial size of predator and prey populations. However, in the stochastic formulation there exists a "catastrophic" sequence of events which leads to depletion of preys by predators and, in consequence, to the extinction of predators as well. When running the model long enough, the probability of not executing this catastrophic sequence drops to zero. This leads to radical qualitative differences in the trajectories obtained by these

A c c e p t e d m a n u s c r i p t

two approaches: in the deterministic case the trajectory in the predator versus prey phase space is an ellipse, while in the stochastic case the trajectory eventually reaches the trivial steady state of no predators and no prey individuals in the system. The expected time it takes to reach this state depends on the initial number of species. Such discrepancy in the trajectories is especially easily observed when the initial population sizes are small.

Another significant impact of random fluctuations can be observed in the model of T cell receptor signalling presented in [START_REF] Lipniacki | Stochastic effects and bistability in T cell receptor signaling[END_REF], where it is shown that, because of bistability of the system and the fact that the T cell activation is due to a small number of foreign peptides, the responses are highly stochastic. This results in stochastic trajectories not following the deterministic trajectory, which converges to a steady state. Instead, the stochastic realisations may occasionally jump between the basins of attraction of two possible states. In particular, as was shown in [START_REF] Lipniacki | Stochastic effects and bistability in T cell receptor signaling[END_REF], stochastic noise can cause a transition from the higher stable state to the lower one and most of the stochastic trajectories are trapped in the basin of attraction of the latter steady state in contrast to the deterministic case.

As a result, the qualitative behaviour revealed by the stochastic approach differs significantly from the behaviour obtained from the deterministic description. For details we refer the reader to [START_REF] Lipniacki | Stochastic effects and bistability in T cell receptor signaling[END_REF].

Although for a complex system detailed mathematical analysis based on the "chemical master equation" is intractable [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF]), it is possible to gain insight into the system's dynamics by performing a series of stochastic simulations of the time-evolution of such system by so called Gillespie algorithm [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF]). The algorithm is a well-established

A c c e p t e d m a n u s c r i p t

procedure for generating a stochastic realisation of the system's temporal behaviour. However, due to reasons such as computational efficiency, availability of dedicated simulation software with analysis tools (steady-state, sensitivity analysis, etc.), and expertise in the theory of differential equations, the deterministic modelling approach is commonly used in examination of biological systems, although the stochastic formulation in many cases would be more justified.

Bearing in mind the above mentioned merits of the new simplified heat shock response model described in Petre et al. (2009a), the aim of this paper is to show that in this particular case approximating a discrete system with a continuous model is a valid approach. A stochastic model complementary to the deterministic one is developed. An effort to perform 1000 stochastic simulations is made in order to investigate whether the qualitative results of the stochastic model agree with the deterministic outcome. Having the problem of small number of molecules of some of the reactants in mind (initial number concentrations of hsf, hsf 2 , hsf 3 , hsf 3 : hse, hse, hsp : mfp, see Petre et al. (2009a) for details), as explained above, one could expect the time-course trajectories obtained with the stochastic model to be substantially different from the trajectories computed in the deterministic formalism.

However, we show that the influence of the random fluctuations does not invalidate the continuous approach and the obtained results support the use of the deterministic formulation in this case. In particular, we investigate the number of steady states of the deterministic model and compare the obtained results with the dynamics demonstrated by the stochastic model. We show that the underlying stochastic process of our model has a unique stationary

A c c e p t e d m a n u s c r i p t

distribution and that the performed stochastic simulation results reveal no evidence of multistationarity, which is consistent with the deterministic description. Additionally, this analysis let us gain some more insight into the dynamics of the heat shock response mechanism. The question about the stationarity and stability, i.e. the number of steady-states and whether they are stable or unstable, is important in the examination of the dynamics of biological systems. For example, bistability in biological systems is, in general, accompanied by hysteresis, which in turn promotes robustness [START_REF] Karmakar | Posititve feedback, stochasticity and genetic competence[END_REF]; [START_REF] Lipniacki | Stochastic effects and bistability in T cell receptor signaling[END_REF]).

The paper is organised as follows. In Section 2 we briefly describe the simplified deterministic model (named deterministic model for compactness in the continuation) of the heat shock response in eukaryotic cells which was proposed in Petre et al. (2009a). Next, in Section 3, we discuss the Markov jump process which constitutes the corresponding stochastic model and show that it has a unique stationary distribution. Further, in Section 4, the stochastic simulation results are discussed and a comparison between the deterministic and stochastic model is presented. Finally, we end with conclusions in Section 5.

Deterministic model

The model of the eukaryotic heat shock response consists of four main modules: the heat-induced protein misfolding, the dynamic transactivation of the genes encoding heat shock proteins, their backregulation and the chaperone activity of the heat shock proteins.

At elevated temperatures proteins tend to misfold and create aggregates,

A c c e p t e d m a n u s c r i p t

which has disastrous effects on the cell. In order to survive, the cell has to promptly increase the level of heat shock proteins (hsp), which is the main task of the heat shock response mechanism. Heat shock proteins act as chaperones: they interact with the misfolded proteins (mfp) and assist them in refolding to their native conformation (prot). The control over the defence mechanism against the temperature-induced harmful phenomena is implemented through the regulation of the transactivation of the hsp-encoding gene. Activation of the transcription proceeds along the following scheme:

heat shock factors (hsf) trimerize (through a transient dimerization) and in this form bind to the heat shock element (hse), i.e. the promoter of the hspencoding gene. Once the hsf-trimer (hsf 3 ) is bound to the specific DNA sequence (hsf 3 : hse), the gene is transactivated and new hsp molecules are eventually synthesised. Finally, when the level of hsps is high enough to cope with the thermal stress, the production is switched off: hsps bind both to free hsfs and hsfs that occur in compound forms (hsf 2 , hsf 3 , hsf 3 : hse), which, in consequence, get disassembled. As a result, DNA transcription of hspencoding gene is turned off and the formation of new hsf trimers is blocked.

The full list of molecular reactions constituting the model is presented in Table 1. By assuming the law of mass-action for all reactions (R 1 )-(R 17 ), the associated mathematical model based on ordinary differential equations is obtained. The rate coefficient of protein misfolding with respect to the temperature (ϕ(T )) in reaction (R 14 ) is given by the following formula:

ϕ(T ) = (1 - 0.4 e T -37 ) • 1.4 T -37 • 1.45 • 10 -5 s -1 , ( 1
)
where T is the numerical value of the temperature of the environment in Celsius degrees. The formula is valid for 37 ≤ T ≤ 45. It is based on

A c c e p t e d m a n u s c r i p t

experimental investigations in [START_REF] Lepock | Protein denaturation in intact hepatocytes and isolated cellular organelles during heat shock[END_REF], [START_REF] Lepock | Thermal analysis of chl v79 cells using differential scanning calorimetry: Implications for hyperthermic cell killing and the heat shock response[END_REF] and was originally proposed in [START_REF] Peper | A mathematical model of the hsp70 regulation in the cell[END_REF]. Expression (1) in its current form was obtained by adapting the original formula to the time unit of the discussed mathematical model (see Petre et al. (2009a)). In our survey the temperature is set to 42 • C, i.e. the cells are exposed to heat shock conditions.

As shown in Petre et al. (2009a), there are three mass-conservation relations in the model: the total number of heat shock factor molecules, heat shock elements and protein molecules (either misfolded or in native conformation) is conserved in time. This can be written formally as

C 1 = hse(t) + 3 hsf 3 : hse(t) (2)
C 2 = hsf(t) + 2 hsf 2 (t) + 3 hsf 3 (t) +3 hsf 3 : hse(t) + hsp : hsf (3) C 3 = prot(t) + mfp(t) + hsp : mfp(t) (4)
for all t ≥ 0, where C 1 , C 2 , C 3 ≥ 0 are some constants determined by initial conditions, i.e. right-hand side expressions at t = 0 in the above equations

(2)-( 4).

The described model of eukaryotic heat shock response is based solely on well-documented reactions and does not include any "artificial" elements such as experimentally unsupported components or biochemical reactions.

For a detailed discussion of the model, we refer the reader to Petre et al. (2009a).

A c c e p t e d m a n u s c r i p t

Stochastic model

Stochastic modelling of biochemical networks is today well-established.

The time-evolution of a reaction system can be regarded as a stochastic process (cf. [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF]). In particular, the dynamics of a biochemical network can be viewed as a continuous-time Markov process. A continuoustime stochastic process {X(t), t ≥ 0} with discrete state space S is said to be a continuous-time Markov chain (CTMC for short) if

P {X(t n ) = i n | X(t 0 ) = i 0 , . . . , X(t n-1 ) = i n-1 } = P {X(t n ) = i n | X(t n-1 ) = i n-1 } for all 0 ≤ t 0 < • • • < t n-1 < t n and i 0 , . . . , i n-1 , i n ∈ S.
The Markov property expresses that the conditional distribution of a future state given the present and past states depends only on the present state and is independent of the past.

We consider a time-homogeneous Markov chain for which the transition

probability P {X(t+u) = j | X(u) = i} is independent of u. Let Q = (q ij ) i,j∈S
be the infinitesimal transition rate matrix of the continuous-time Markov chain {X(t)} such that the following assumption is satisfied.

Assumption 1. ν i = j =i q ij are positive and bounded in i ∈ S.

As stated in [START_REF] Tijms | Chaperone regulation of the heat shock protein response[END_REF]

p ij = ⎧ ⎨ ⎩ q ij /ν i j = i, 0 j = i,
for all i, j ∈ S (e.g. see [START_REF] Tijms | Chaperone regulation of the heat shock protein response[END_REF]).

The corresponding stochastic formulation of the HSR molecular model presented in Section 2 is a system of 10 chemically active species S i , (i = 1, . . . , 10), (S 1 ≡ hse, S 2 ≡ hsf, S 3 ≡ hsf 2 , S 4 ≡ hsf 3 , S 5 ≡ hsf 3 : hse, S 6 ≡ hsp, S 7 ≡ hsp : hsf, S 8 ≡ hsp : mfp, S 9 ≡ mfp, S 10 ≡ prot) that participate in 17 chemical reactions (R 1 ) -(R 17 ) in some volume V . More specifically, in our case the abstract volume V is simply a eukaryotic cell.

A c c e p t e d m a n u s c r i p t

The system's state space S ⊂ N 10 is defined as

S = {(s 1 , . . . , s 10 ) T : ∀ i∈{1,...,10} s i ∈ N, s 1 + s 5 = C 1 , s 2 + 2 s 3 + 3 s 4 + 3 s 5 + s 7 = C 2 , s 8 + s 9 + s 10 = C 3 , s 6 + s 7 + s 8 ≤ K}, ( 5
)
where C 1 ≥ 1, C 2 ≥ 3 and C 3 ≥ 1 are constants describing the fixed in time total number of hse, hsf and protein (either misfolded or in native conformation) molecules present in the system, respectively. The last inequality, i.e. s 6 + s 7 + s 8 ≤ K, requires some comment. In the absence of it the model would make allowance for any unbounded number of free hsp molecules to co-exist. However, this is certainly contrary to the fact that any living cell has a limited volume and, in consequence, can only contain a finite number of hsp molecules. Thus, in order to make the model more realistic, a big enough (in the sense that it allows for the appropriate number of free hsp molecules to be present in the system) constant K is introduced and an upper bound on the value of the S 6 variable is imposed by the last inequality. The direct consequence of adding it is that the state space S becomes finite.

The system is in state s = (s 1 , . . . , s 10 ) T at time t if and only if the number of molecules of species S i at time t is s i for all i = 1, . . . , 10. The conditions posed on the constants C 1 , C 2 and C 3 ensure that the HSR mechanism is operational, i.e. that at least one hse molecule is present in the system, that the system is able to produce at least one hsf 3 molecule that can bind to the DNA and, in consequence, initiate the transcription and translation of hsp.

Finally, that at least one generic protein prone to misfolding exists in the Each reaction R μ is characterised by a stochastic rate constant c μ (μ = 1, . . . , 17), see Table 3. The values were obtained from the deterministic model (Petre et al., 2009a). In the case of reaction R 1 , the deterministic rate constant value was multiplied by 2 in order to obtain the value for the corresponding stochastic rate constant.

Let v μ be the μ-th column of the stoichiometry matrix of the HSR system presented in Table 2. Reaction R μ causes the system to make a transition from some state i ∈ S to state j = i + v μ . The fundamental hypothesis of the stochastic formulation of chemical kinetics [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF] is that the reaction parameter c μ can be defined as follows.

c μ dt ≡ average probability, to first order in dt, that a particular combination of R μ reactant molecules will react accordingly in the next time interval dt.

As shown by Gillespie in [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], c μ is dependent on the radii of the molecules involved in the reaction and their average relative velocities, where the average relative velocity is a function of the temperature of the system and the individual molecular masses. Further, it is shown that the probability of reaction R μ occurring in V in the time interval (t, t +dt), given that the system is in state i at time t, has the form h i μ c μ dt. h i μ denotes the number of possible combinations of reactant molecules involved in reaction R μ when the system is in state i. However, since the total number of hsp molecules that might co-exist in a cell is limited, no further hsp production (reaction (R 7 )) should take place when the system is in any of the states in A c c e p t e d m a n u s c r i p t which the limit is reached. Let us denote all these states by S K , i.e.

S K = {s | s 6 + s 7 + s 8 = K}. (6)
Thus, the probability of reaction R 7 occurring in V in the time interval (t, t + dt) when the system is in i ∈ S K should be 0. Hence

h i 7 = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ No. of combinations of reactant molecules of R 7 i ∈ S \ S K , 0 i ∈ S K . ( 7
)
Due to the fact that the reaction hazards depend only on the current state of the system, the time-evolution of the state of the reaction system can be regarded as a CTMC. Since the state space S is finite, Assumption 1 is fulfilled and hence the chain is a Markov jump process constructed as described above. The infinitesimal transition rates of the Markov jump process are

q ij = h i μ c μ , (8)
where

j = i + v μ .
Now, let us consider the reaction probability density function P i (τ, μ) dτ of the HSR system, which forms the basis for the Gillespie's simulation algorithm. It is defined by

P i (τ, μ) dτ ≡ probability at time t (when the system is in state i ∈ S)
that the next reaction in V will occur in the differential time interval (t + τ, t + τ + dτ ), and will be an R μ reaction.

A c c e p t e d m a n u s c r i p t

As shown in [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], the density function can be expressed as

P i (τ, μ) dτ = h i μ c μ • exp - 17 ν=1 h i ν c ν τ dτ . ( 9
)
The superscript i ∈ S indicates that in fact we deal with a whole family of such functions. Which of them is considered at time t depends on the state of the system at time t. In the continuation, in order to lighten the language, "probability at time t" will be a shorthand for "probability at time t when the system is in state i". Let P i 1 (τ ) dτ denote the probability at time t that the next reaction will occur between times t + τ and t + τ + dτ , irrespective of which reaction it might be. By the definition of P i (τ, μ) we have that

P i 1 (τ ) = 17 μ=1 P i (τ, μ) (10) = ( 17
k=1 h i k c k ) • exp[-( 17
k=1 h i k c k ) • τ ].
Hence, the sojourn-time rates ν i of the Markov jump process are given by

ν i = 17 k=1 h i k c k . ( 11
)
The probability of the transition from state i to state j = i+v μ of the Markov jump process (and, in consequence, of the embedded Markov chain {X n }) is the probability at time t that the next reaction in V will be an R μ reaction.

Using equation ( 9), it can be expressed as

p ij = ∞ 0 P i (τ, μ) dτ = h i μ c μ 17 k=1 h i k c k (12)
if j = i and p ii = 0 for all i ∈ S.

Lemma 1. The embedded Markov chain {X n } is irreducible.

A c c e p t e d m a n u s c r i p t

Before presenting the proof let us divide the species into two groups. The first one, called elementary species group (denoted by G elementary ), contains hsf, hse, hsp, mfp and prot species. The other one, denoted by G compound and named compound species group, is made of all the remaining species.

Proof of Lemma 1. Let i, j be any two states from the state space S. By i j we denote that the state j is reachable from the state i, i.e. that there exists a sequence of reactions (R 1 )-(R 17 ) which leads the system from the state i to the state j.

In order to prove that {X n } is irreducible it is enough to show that i j, since i and j are two arbitrarily chosen states. Let hsp(k), mfp(k) and prot(k)

be the total number of hsp, mfp and prot molecules present in the system when in the state k, respectively. Further, let z = (C 1 , C 2 , 0, 0, 0, hsp(i), 0, 0, mfp(i), prot(i)) T .

z is obtained from i by disassembling all compound species from G compound .

Thus, in the state z the number of molecules of any species from G compound is 0 and the number of molecules of any s ∈ G elementary is equal to the total number of s molecules in the system in the state i. Clearly z ∈ S and i z, since for any s ∈ G compound there exists a sequence of reactions (R 1 )-(R 17 )

which disassembles s into elements from G elementary .

Let z = (C 1 , C 2 , 0, 0, 0, hsp(j), 0, 0, mfp(i), prot(i)) T and z = (C 1 , C 2 , 0, 0, 0, hsp(j), 0, 0, mfp(j), prot(j)) T .

A c c e p t e d m a n u s c r i p t

z ∈ S and we continue to show that z z . There are three cases: hsp(i) = hsp(j), hsp(j) < hsp(i) or hsp(j) > hsp(i). In the first case z = z and trivially z z . If hsp(j) < hsp(i), z can be reached from z by applying reaction (R 13 ) hsp(i)hsp(j) times. If finally hsp(j) > hsp(i), first hsf 3 : hse is produced (this is doable since C 1 ≥ 1 and C 2 ≥ 3). Next, by applying reaction (R 7 ) hsp(j)hsp(i) times the required number of additional hsp molecules is produced. Finally, hsf 3 : hse is disassembled by applying a sequence of reactions <(R 6 ), (R 3 ), (R 2 )>. Hence z z .

We continue to show that z z . There are two cases. Either mfp(i) > mfp(j) or prot(i) ≥ prot(j) since mfp(k) + prot(k) = C 3 for any state k ∈ S.

In the first case, if hsp(j) = 0, first one hsp molecule is produced by applying reaction sequence <(R 1 ), (R 3 ), (R 5 ), (R 7 ), (R 6 ), (R 4 ), (R 2 )>, which leads to state (C 1 , C 2 , 0, 0, 0, 1, 0, 0, mfp(i), prot(i)) T . Then, by applying reaction sequence <(R 15 ), (R 17 )> mfp(i)mfp(j) times the system reaches state:

(C 1 , C 2 , 0, 0, 0, 1, 0, 0, mfp(j), prot(j)) T since, as mentioned before, mfp(k) + prot(k) = C 3 for any state k ∈ S.

Finally, the only hsp molecule is degraded by applying reaction (R 13 ) and the system arrives in state z . If hsp(j) is greater than 0, state z can be reached by less steps since neither production nor degradation of the one additional hsp molecule is required.

In the second case, when prot(i) ≥ prot(j), state z can be reached by applying misfolding reaction (R 14 ) prot(i)prot(j) times. Hence z z .

At last z j. State j is reached by producing the appropriate numbers of molecules of all compound species. Since in state z the required number of molecules of all elementary species is already present, by apply-

A c c e p t e d m a n u s c r i p t

ing appropriate reactions all compound species molecules can be produced.

The mass-conservation law ensures that the numbers of elementary species molecules will be decreased appropriately and that the system will reach state j. Hence {X n } is irreducible.

The irreducibility of the embedded chain {X n } implies the irreducibility of the continuous-time Markov chain {X(t)}. Since the state space S is finite, it follows that the CTMC {X(t)} is positive recurrent. In consequence, it has an invariant measure η which is unique up to multiplicative factors and can be found as the solution of the equation

η T Q = 0. Moreover, i∈S η i < ∞
since S is finite and there exists a unique stationary distribution π of {X(t)}

given by

π = η i k∈S η k i∈S . ( 13
)
For the theoretical details we refer the reader to, e.g., [START_REF] Norris | Markov Chains. Cambridge Series in Statistical and Probabilistic Mathematics[END_REF]Resnick (1992).

Results and discussion

The deterministic approach, based on the law of mass action, yields a system of ordinary differential equations for molecular concentrations. In consequence, the biochemical system is modelled as being continuous. But such description does not capture effects that occur due to either the discreteness of molecular quantities or the stochastic nature of chemical reactions of some reactants become smaller [START_REF] Mcadams | It's a noisy business! Genetic regulation at the nanomolar scale[END_REF][START_REF] Srivastava | Stochastic vs. deterministic modeling of intracellular viral kinetics[END_REF]. This is the case of the deterministic heat shock response model being discussed: except for prot, hsp, hsp : hsf and mfp, all the other species have very small initial number of molecules (Table 3) and, as can be seen from the continuous simulation results, stay at the low level throughout the time of simulation. This might be the main objection to the continuous approach applied in Petre et al. (2009a,b). Since the stochastic modelling seems more reasonable in this case, we made the effort to run 1000 stochastic simulations in order to check whether the dynamics of the continuous description agrees qualitatively with the behaviour demonstrated by the discrete system. The results of 1000 independent stochastic simulation runs (blue and green points) for 5 species: hsf 3 : hse, hsp, mfp, hsp : hsf and hsp : mfp, overlaid with the deterministic outcome (yellow line) are shown in Figure 1.

The mean together with the mean +/-standard deviation are shown in Figure 2. The ratios of the sample standard deviation to the sample mean were computed for the 5 considered species and are depicted in Figure 3.

According to [START_REF] Gillespie | A general method for numerically simulating the stochastic time evolution of coupled chemical reactions[END_REF], since the ratios are small (less than 0.12 in the case of mfp and hsp, see Fig. 3c and3b) and very small (less than 0.035 for hsp : mfp and less than 0.007 for hsp : hsf, see Fig. 3e and3d, respectively), the results of independent runs of the system are expected not to vary much and the presented outcomes of 1000 stochastic simulations together with the estimated mean should provide a statistically adequate picture of the evolution of the chemical system in time. One might argue that the ratio for hsf 3 : hse is however quite big: it peaks at about 1.1 and stabilises below 0.6 (see Fig. 3a). In this particular case the mean converges to approximately 3

A c c e p t e d m a n u s c r i p t

molecules of hsf 3 : hse and the standard deviation is around 1.6, which all in all gives a narrow range of possible values of molecule number and hence this result can be accepted.

We first investigated the number of steady-states of the deterministic model. Since our attempts to analytically solve the algebraic system of steady state equations obtained from the differential ones did not bring any results, we performed some numerical investigations. We randomly chose 10000 sets of initial particle numbers for the continuous model from a wide range of values, but in such a way that the total amounts of hse, hsf and proteins in the resulting system would always be the same as in the case of the original deterministic model presented in Petre et al. (2009a). For each of these sets we run numerical time-course simulations and waited for the considered system to stabilize. In all these cases the systems converged to exactly the same state as the original model, i.e. no other steady states were found by this method.

Additionally, bifurcation analysis performed with the AUTO software (XP-PAUT was used as the front-end, [START_REF] Doedel | AUTO 97, software for continuation and bifurcation problems in ordinary differential equations[END_REF]; [START_REF] Ermentrout | Simulating, Analyzing, and Animating Dynamical Systems: A Guide to XPPAUT for Researchers and Students[END_REF])

with respect to parameter values did not reveal multistationarity (data not shown). These results suggest that the heat shock response mechanism is rather monostable.

Next, we were interested in investigating the range of behaviour the stochastic model was likely to exhibit. As shown in Section 3, there exists only one stationary limit distribution π given by Eq. ( 13), which governs the transitions of the Markov jump process when the number of iterations goes to infinity. In particular, we analysed the unimodality of the hsp level by computing some appropriate statistics from the performed 1000 stochastic Next, in order to check whether the realisations r i , i = 1, . . . , 1000, can be divided into subgroups such that the means of the subgroups would differ significantly from each other, we applied the following procedure. We defined two subsets:

S U = {r i : ∀ t∈T r i (t) > m(t) - s 4 ∧ ∃ t∈T r i (t) > m(t) + s 4 and S L = {r i : ∀ t∈T r i (t) < m(t) + s 4 ∧ ∃ t∈T r i (t) < m(t) - s 4 .
In Further, a clustering algorithm was applied in order to determine whether some subsets of realisations could be isolated and the computed means would point to potential multimodality. To this aim, we utilised the Agnes algorithm (implementation of an agglomerative hierarchical clustering method, [START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF]) with the manhattan metric, i.e. the distance between two realisations r i and r j is defined as d(r i , r j ) = t∈T |r i (t) - Finally, as suggested in [START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF], we investigated the empirical probability mass function by drawing histograms of the realisations at some time point in the considered time interval T . Figure 6 shows the histograms overlaid with the normal distribution curve with mean and standard deviation computed from all 1000 realisations at time point t = 150000s. In the case of Fig. 6a, where the bin width is set to 20, the obtained results indicate that the distribution is unimodal. Changing the bin size to 10 (Fig. 6b) does not change the picture significantly.

Although due to the small particle numbers of some of the reagents the stochastic modelling is more reasonable, the presented results do not reveal

Conclusions and Further research

In this paper we presented a stochastic model associated with a previously described (Petre et al., 2009a) model of the heat shock response in eukaryotic cells. The stochastic model was viewed as a Markov jump process and the existence and uniqueness of the stationary distribution was shown.

Further, the model was compared to the deterministic description of heat shock response (Petre et al., 2009a). The aim with the comparison was to show that in this particular case the approximation of a discrete system with a continuous model is reasonable. This is not true in general, especially when the numbers of metabolites in the considered biochemical system are small.

The presented results indicate that the stochastic and deterministic models provide a qualitatively consistent picture of the dynamics of the heat shock response mechanism. Additionally, the development of the stochastic model and the effort of performing 1000 stochastic simulations enabled gaining some more information about the dynamics of the heat shock response. The outcomes of the analysis of the stochastic realisations lead towards the conclusion that the heat shock response mechanism is a rather monostable system.

Moreover, this is in agreement with the results of the analysis performed on the deterministic model. All in all, the presented results strengthen the trust in the deterministic description of the HSR mechanism in eukaryotic cells

A c c e p t e d m a n u s c r i p t

proposed in Petre et al. (2009a).

Although it was shown in Section 3 that the Markov jump process has a unique stationary distribution, there is no certainty that it was reached already in the considered time interval T = {130000s, . . . , 150000s}. It was chosen based on the results of many stochastic simulations, which suggest that the process stabilises relatively long before the time point t = 130000s.

Nevertheless, some assessment of the convergence to the stationary distribution in this case would be desired. One of possible approaches is to measure the rate of convergence by the mixing time [START_REF] Sinclair | Improved bounds for mixing rates of markov chains and multicommodity flow[END_REF]. For ergodic Markov chains the rate is governed by the second largest eigenvalue in absolute value λ 2 , in particular the spectral gap 1λ 2 is both a necessary and sufficient condition for rapid mixing, see [START_REF] Sinclair | Improved bounds for mixing rates of markov chains and multicommodity flow[END_REF]

2 hsf → hsf 2 (R 1 ) hsf 2 → 2 hsf (R 2 ) hsf + hsf 2 → hsf 3 (R 3 ) hsf 3 → hsf + hsf 2 (R 4 ) hsf 3 + hse → hsf 3 : hse (R
⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

  other formalisms: we develop a stochastic model associated with the simplified version of the model fromPetre et al. (2009b) which has been described

(

  [START_REF] Mcadams | It's a noisy business! Genetic regulation at the nanomolar scale[END_REF];[START_REF] Pahle | Biochemical simulations: stochastic, approximate stochastic and hybrid approaches[END_REF];[START_REF] Sandmann | Discrete-time stochastic modeling and simulation of biochemical networks[END_REF][START_REF] Wilkinson | Stochastic Modelling for Systems Biology[END_REF]). As discussed in Section 1, random fluctuations may have a significant impact on the reaction dynamics, especially as the numbers of molecules

  computed the median m(t) of the 1000 stochastic realisations on the time interval T = {130000s, . . . , 150000s}. It is depicted in Fig. 4 as the middle black line. The upper and lower black lines are m(t) +/-1 4 • s respectively, where s is the range of dynamics the model exhibits in the 1000 realisations on the considered time interval, i.e. s = max t∈T, i∈I {r i (t)}min t∈T, i∈I {r i (t)}, where I = {1, . . . , 1000} and r i is the i-th realisation. The mean (brown line) basically coincides with the median on the whole time interval.

  our case, there are 253 realisations in S U and 189 in S L . The means computed from the realisations of each of these subsets are depicted by red lines in Fig.4. The means are close to the global mean on the whole considered time interval and since the numbers of elements in the S U and S L subsets , i.e. approximately 1/4 and 1/5 of all the 1000 considered realisations, this result does not indicate any significant split.

  r j (t)|, thus the realisations are treated as points in a |T |-dimensional space. By applying this metric the characteristics of the realisations on the whole considered time interval are taken into account, hence they are compared in a "global" sense. The obtained dendrogram is presented in Fig. 5. The agglomerative coefficient (AC), which measures the clustering structure of the dataset, is 0.82. This indicates that the clustering algorithm did find some rather clear structuring 1 . We isolated two groups of realisations that stand out on the obtained dendrogram. They are marked in Fig. 5 by two rectangles which enclose the dendrogram branches constituting these groups. The two resulting subclusters are at almost the same height in the clustering tree. The means of the stochastic realisations belonging to these two groups at time point t = 150000s are 757 (left subcluster) and 794 (right subcluster).Although the agglomerative coefficient indicates some clustering structure of the realisations, the mean values are very close to each other and agree well state value of the deterministic model (767).

  any qualitative discrepancy in the dynamics of the two considered models of the heat shock response. The range of behaviour the stochastic model is likely to exhibit, which can be observed based on the performed 1000 simulations, confirms the dynamics of the continuous model. The performed analysis of the stochastic realisations does not reveal any clear signs of multistationarity of the HSR mechanism. Although unimodality of a stationary probability density function does not necessarily imply the uniqueness of the stable steady state of the deterministic approximation (as well as bimodality does not determine the existence of bistability, etc.), usually this is the case and to this extent the stochastic results agree with the deterministic outcomes indicating that there exists only one stable steady state. This shows that the approximation of a discrete system with a continuous model is valid and strengthens the trust in the deterministic description. Additionally, the presented stochastic formulation, together with the performed analysis of its to the continuous description, let us gain more insight into the dynamics of the HSR mechanism, especially in respect of the number of steady states, which, as discussed previously, is important from a biological point of view.

R

  Development Core Team, 2008. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Aus-

Figure 3 :

 3 Figure 3: The ratios of the standard deviation to the sample mean at each considered time point.

Figure 4 :

 4 Figure 4: The median of the 1000 realisations on the time interval T = {130000s, . . . , 150000s} (middle black line). The upper and lower black lines are the median ±1/4 of the range of dynamics the model exhibits in the 1000 realisations on the considered time interval. The mean of all the realisations, of the subset S U and S L plotted with brown, upper red and lower red lines, respectively.

  ) if the system jumps to state i, it then stays in state i for an exponentially distributed time with mean 1/ν i independently of how the system reached state i and how long it took to get there (this explains the name sojourn-time rates used for {ν i }); b) if the system leaves state i, it jumps to state j (j = i) with probability p ij independently of the duration of the stay in state i.

	, if Assumption 1 is fulfilled, it can be shown A c c e p t e d m a n u s c r i p t
	that the infinitesimal transition rates determine a unique continuous-time
	Markov chain which is precisely a Markov jump process constructed in the
	following way:

aFurther, let {X n , n = 0, 1, . . . } be the embedded Markov chain, i.e. X n is defined as the state of {X(t)} just after the n-th transition with the convention that X 0 = X(0). The one-step transition probabilities of the discrete-time Markov chain {X n } are given by

  for details. The problem of determining λ 2 of the presented Markov chain underlying the stochastic model of heat shock response is subject of further research.

	imental biology. (Eds.), Mathematical Theory and Computational Practice: 5th Confer-
	ence on Computability in Europe, CiE 2009, Heidelberg, Germany, July
	19-24, 2009, Proceedings. Vol. 5635 of Lecture Notes in Computer Science.
	Springer, pp. 399-408.
	The rate constant values for the presented stochastic model were ob-tained from the corresponding values of the deterministic model presented in Petre et al. (2009b), which in turn were fitted to available experimen-tal data. As suggested in Wilkinson (2006), another way of deducing the rate constant values for the stochastic model could utilise methods that are based on Bayesian inference and take advantage of Markov Chain Monte Carlo (MCMC) algorithms such as the Metropolis-Hastings algorithm or the Gibbs Sampler. However, such methods demand high-quality, calibrated, A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t Petre, I., Mizera, A., Hyder, C. L., Mikhailov, A., Eriksson, J. E., Sistonen, L., Back, R.-J., 2009b. A new mathematical model for the heat shock response. In: Condon, A., Harel, D., Kok, J. N., Salomaa, A., Winfree, E. (Eds.), Algorithmic Bioprocesses. Natural Computing Series. Springer, pp. 411-425. Pockley, A., 2003. Heat shock proteins as regulators of the immune response. The Lancet 362 (9382), 469-476. a n u s c r i p t Powers, M., Workman, P., 2007. Inhibitors of the heat shock response: Biol-m ogy and pharmacology. FEBS Lett. 581 (19), 3758-3769. Remondini, D., Bernardini, C., Forni, M., Bersani, F., Castellani, G., Bacci, M., 2006. Induced metastable memory in heat shock response. J. Biol. Physics 32, 49-59. Resnick, S., 1992. Adventures in Stochastic Processes. Birkhäuser. Rieger, T. R., Morimoto, R. I., Hatzimanikatis, V., 2005. Mathematical mod-A c c e p t e d eling of the eukaryotic heat shock response: Dynamics of the hsp70 pro-
	high-resolution time-course measurements for a reasonably large subset of moter. Biophysical Journal 88 (3), 1646-58.
	model metabolites (Wilkinson, 2006). Unfortunately, experimental data of
	such quality are still seldom if ever available and make a challenge for exper-

Table 1 :

 1 The simplified model for the eukaryotic heat shock response

	5 )

Table 2 :

 2 The stoichiometry matrix of the heat shock response model. The columns correspond to reactions (R 1 )-(R 17 ) and the rows to metabolites in the order: hse, hsf, hsf 2 , hsf 3 , hsf 3 : hse, hsp, hsp : hsf, hsp : mfp, mfp, prot.

	Param. Reaction Value	Unit	Metabolite Init. no.
	k + 1	(R 1 )	6.98	V #•s	hsf	0
	k -1	(R 2 )	0.19	s -1	hsf 2	0
	k + 2	(R 3 )	1.07	V #•s	hsf 3	0
	k -2	(R 4 )	10 -9	s -1	hse	29
	k + 3 k -3 k 4 k + 5 k -5 k 6 k 7 k 8 k 9 k 10 k + 11 k -11 k 12	(R 5 ) (R 6 ) (R 7 ) (R 8 ) (R 9 ) (R 10 ) (R 11 ) (R 12 ) (R 13 ) (R 14 ) (R 15 ) (R 16 ) A c c e p t e d 0.17 1.21 • 10 -6 8.3 • 10 -3 9.74 3.56 2.33 4.31 • 10 -5 2.73 • 10 -7 3.2 • 10 -5 ϕ(42) = 7.77 • 10 -5 s -1 V #•s s -1 s -1 V #•s s -1 V #•s V #•s a n u s c r i p t hsf 3 : hse hsp 766 2 hsp : hsf 1403 mfp 517 hsp : mfp 71 prot 1.15 • 10 8 V #•s m s -1 3.32 • 10 -3 V #•s 4.44 s -1 (R 17 ) 13.94 s -1

Table 3 :

 3 The numerical values of the parameters and the initial numbers of molecules in the stochastic model. The numerical quantities are obtained by adopting the corresponding values in Petre et al. (2009a): the initial numbers of molecules are truncated to natural numbers, the value of the rate constant k + 1 is twice the value of the corresponding deterministic rate constant. # denotes the number of molecules, V is the cell volume and s -second.

AC is a dimensionless quantity, varying between 0 and 1 -AC close to 1 shows that a very clear structure has been found, while value 0 implies that the data consists of only one big cluster, see e.g.[START_REF] Kaufman | Finding Groups in Data: An Introduction to Cluster Analysis[END_REF] for details.

Acknowledgements

Both deterministic and stochastic models were implemented and run in Copasi, a software application for simulation and analysis of biochemical networks (Hoops et al., 2006). The stochastic simulations were performed using the Gibson and Bruck algorithm (Gibson and Bruck, 1998). The obtained time-course data were analysed and plotted in R, a software environment for statistical computing and graphics (R Development Core Team, 2008).

The authors give special thanks to Anna Gambin from the Faculty of Mathematics, Informatics and Mechanics, University of Warsaw for helpful advice as well as thorough and detailed comments. Andrzej Mizera would like to express his gratefulness to Jan Westerholm, Mats Aspnäs and Evren Yurtesen from the Department of Information Technologies, Åbo Akademi University for making available their computational resources and providing essential help for performing the stochastic simulations. The author would like to thank Ion Petre from the Computational Biomodelling Laboratory, Åbo Akademi University for invaluable discussions.

The authors were partially supported by Ministry of Science and Higher Education (grant number N N518 426936).


Romain Delpoux
email: romain.delpoux@ec-lille.fr

Thierry Floquet
email: thierry.floquet@ec-lille.fr

On-line parameter estimation of a magnetic bearing

Keywords: Parameter estimation, Laplace transform, Magnetic bearing. I

This article presents a parameter estimation algorithm for a magnetic bearing. Such process are inherently unstable systems with strongly nonlinear dynamics. Here, a simplified model of the magnetic bearing is developed. This enables to obtain a linear expression with respect to the unknown parameters. These parameters are measurable with difficulties, and may slightly vary over time. The expression of the estimates is written as a function of integrals of the inputs and outputs of the system. The simulations and the experiments show a fast and robust on-line identification.

when identifying parameters. In Section III, expressions for the on-line identification of the desired parameters are developed. Here is treated the on-line estimation of parameters depending on the radial bearing geometry and the materials. The last Section is devoted to simulation and experimental results.

II. PROBLEM STATEMENT

In this article is presented an approach to estimate parameters of the radial bearing (the axial bearing case will be considered in future work). The parameters that are to be estimated here depend on the geometry and the material of the bearing. The estimation of these parameters is important because such parameters are difficult to calculate and may slightly vary over time. The rotor is levitating using a three-phase electromagnetic radial bearing, arranged like three coupled "horseshoe magnets" around the rotor (Fig. 1). The three generated current provide three independent control inputs. The mathematical model of the bearing is based on the assumption of a rigid body and leads to decoupled equations for forces. The dynamics equations, under simplifying assumptions, are written as follows:

m Ÿ = F y (1) m Z = F z (2)
where Y and Z represent the coordinates of the center of mass of the rotor in a Cartesian frame (with axes y et z) which is fixed in the space, at a point being considered as the center of the device. The forces F y and F z represent the resulting forces applied in directions y et z, respectively. The rotor has a mass m. The resultant forces in the plan (yz) are given by the superposition of the forces generated by the Fig. 1. Sketch of the three-phase radial bearing magnets:

F y F z = sin α 1 sin α 2 sin α 3 cos α 1 cos α 2 cos α 3   F 1 F 2 F 3   (3)
The angles that appear in (3) are presented Fig. 1. Individually, the magnetic forces can be modeled by (k ∈ {1, 2, 3})

F k = λ k i 2 k s - sin α k cos α k T Y b Z b 2 (4)
where Y b and Z b are the positions in the bearing plan. s is the nominal air gap and λ k are parameters depending on the geometry and the materials of the bearing. They will be estimated on-line.

Calculating control currents:

The reference currents are obtained from the desired forces. In order to simplify the notations, we assume that the bearing is symmetric, i.e. α 1 = π, α 2 = -π 3 and α 3 = π 3 . In each of the models detailed above, we obtain a couple of independent forces for each bearing plan. Define F 1 such that:

F 1 = F 0 if F z ≥ |Fy| √ 3 F 0 - |Fy| √ 3 + F z else (5)
With F 0 ≥ 0 arbitrarily chosen, the two remaining forces are obtained from the model (3):

F 2 = F 1 + F y √ 3 + F z F 3 = F 1 - F y √ 3 + F z
Then the currents are calculated from the desired magnetic forces using relation [START_REF] Eckhardt | High precision synchronous tool path tracking with an amb machine tool spindle[END_REF]. The above equations reflect the dynamics of the unperturbed model. Model simplifications in (1)-( 2) and modeling errors (e.g. due to incorrect relationship between the input currents i * and the corresponding forces F * ) are assumed to have the same effect than perturbations. The rotor rotation is also a perturbation source with harmonic components. When the rotor moves along an ellipsoid trajectory, the positions in the directions y et z are sinusoidal functions. A significant error in the relation between the currents and the forces can also change sinusoidally. Therefore, the perturbations must be rejected in the λ i estimation algorithm (or must be estimated in order to be compensated in the control).

The perturbed model can be written as:

m Ÿ = F y + p y (t) (6) m Z = F z + p z (t) (7)
where p * (t) represent the perturbations on each axis. Using equation ( 6), ( 7), ( 3) and ( 4), we obtain a relationship of the form:

Ÿ = 3 k=1 λ k u y,k + p y (t) m (8) Z = 3 k=1 λ k u z,k + p z (t) m (9) with u y,k = sin(α k )i 2 k m 0 @ s- 0 @ sin α k cos α k 1 A T 0 @ Y b Z b 1 A 1 A 2 and u z,k = cos(α k )i 2 k m 0 @ s- 0 @ sin α k cos α k 1 A T 0 @ Y b Z b 1 A 1 A 2
These two equations depend on the measured inputs and the accelerations in the plan. The parameters λ k are the parameters to be estimated from measurements Y and Z and despite the perturbations. They may be difficult to compute. The identification of magnetic bearings is usually slow to be implemented and is often made off-line. In the next Section, we will develop expressions for the identification leading to an on-line, fast and robust estimation of these parameters.

III. ALGEBRAIC APPROACH

A. Estimation of λ k parameters

It will be shown in this Section, that the three parameters λ 1 , λ 2 and λ 3 can be estimated using equation [START_REF] Levine | A nonlinear approach to the control of magnetic bearings[END_REF]. Indeed, using an algebraic approach for this system, it is possible to express the parameters λ k as a function of the measured outputs and the inputs only. In order to simplify the computations, it is assumed that the perturbations are constant. We will see in the next paragraph how to deal with non constant perturbations. Consider the Laplace transform of (9) :

s 2 Z(s)-sZ(0)-Ż(0) = λ1UZ,1(s)+λ2UZ,2(s)+λ3UZ,3(s)+ a s (10)
and multiply (10) by s:

s 3 Z(s) -s 2 Z(0) -s Ż(0) = λ 1 sU Z,1 (s) + λ 2 sU Z,2 (s) + λ 3 sU Z,3 (s) + a (11)
Derive three times with respect to s to cancel the initial conditions and the perturbations:

6Z(s) + 18s d ds (Z(s)) + 9s 2 d 2 ds 2 (Z(s)) + s 3 d 3 ds 3 (Z(s)) = λ1 " 3 d 2 ds 2 (UZ,1(s)) + s d 3 ds 3 (UZ,1(s)) " +λ2 " 3 d 2 ds 2 (UZ,2(s)) + s d 3 ds 3 (UZ,2(s)) " +λ3 " 3 d 2 ds 2 (UZ,3(s)) + s d 3 ds 3 (UZ,3(s)) " (12)
We remind that differentiation with respect to s in the operational domain results in a multiplication by -t in the time domain. The multiplication by s in the operational domain leads to the derivation in the time domain. The application of the linear estimator ( 12) is then not convenient. Derivation amplify the high frequency and then the noise contribution. A simple solution is to make the estimator proper. We simply divide ( 12) by s 4 in order to cancel the derivation terms and obtain a relationship with integral operators:

6s -4 Z(s) + 18s -3 d ds (Z(s)) + 9s -2 d 2 ds 2 (Z(s)) + s -1 d 3 ds 3 (Z(s)) = λ1 " 3s -4 d 2 ds 2 (UZ,1(s)) + s -3 d 3 ds 3 (UZ,1(s)) " +λ2 " 3s -4 d 2 ds 2 (UZ,2(s)) + s -3 d 3 ds 3 (UZ,2(s)) " +λ3 " 3s -4 d 2 ds 2 (UZ,3(s)) + s -3 d 3 ds 3 (UZ,3(s)) " (13)
To return to the time domain, compute the inverse Laplace transform and apply the following formula:

Z t 0 Z t k-1 0 ... Z t 1 0 x(τ )dτ dt 1 ...dt k-1 = 1 (k -1)! Z t 0 (t-τ ) k-1 x(τ )dτ (14)
One gets:

R t 0 ˆ6 3! (t -τ ) 3 -18 2 (t -τ ) 2 τ + 9(t -τ )τ 2 -τ 3 ˜Z(τ )dτ = λ1 " R t 0 ˆ3 3! (t -τ ) 3 τ 2 -1 2 (t -τ ) 2 τ 3 ˜uZ,1(τ )dτ " + λ2 " R t 0 ˆ3 3! (t -τ ) 3 τ 2 -1 2 (t -τ ) 2 τ 3 ˜uZ,2(τ )dτ " + λ3 " R t 0 ˆ3 3! (t -τ ) 3 τ 2 -1 2 (t -τ ) 2 τ 3 ˜uZ,3(τ )dτ " (15)
The previous equation depends on the integrals of the measured outputs and inputs. However, there is three unknown parameters for one equation. Then, we have to generate two more equations in order to have the same number of equations than unknowns. To obtain a second equation, take equation ( 11) but differentiating it four times with respect to s and apply the same manipulations as before. We obtain:

R t 0 ˆ-24 3! (t -τ ) 3 τ + 36 2 (t -τ ) 2 τ 2 -12(t -τ )τ 3 + τ 4 ˜Z(τ )dτ = λ1 " R t 0 ˆ-4 3! (t -τ ) 3 τ 3 + 1 2 (t -τ ) 2 τ 4 ˜uZ,1(τ )dτ " +λ2 " R t 0 ˆ-4 3! (t -τ ) 3 τ 3 + 1 2 (t -τ ) 2 τ 4 ˜uZ,2(τ )dτ " +λ3 " R t 0 ˆ-4 3! (t -τ ) 3 τ 3 + 1 2 (t -τ ) 2 τ 4 ˜uZ,3(τ )dτ " (16)
The last equation is obtained, from equation ( 11) but differentiating it five times with respect to s. One gets:

R t 0 ˆ60 3! (t -τ ) 3 τ 2 -60 2 (t -τ ) 2 τ 3 + 15(t -τ )τ 4 -τ 5 ˜Z(τ )dτ = λ1 " R t 0 ˆ5 3! (t -τ ) 3 τ 4 -1 2 (t -τ ) 2 τ 5 ˜uZ,1(τ )dτ " +λ2 " R t 0 ˆ5 3! (t -τ ) 3 τ 4 -1 2 (t -τ ) 2 τ 5 ˜uZ,2(τ )dτ " +λ3 " R t 0 ˆ5 3! (t -τ ) 3 τ 4 -1 2 (t -τ ) 2 τ 5 ˜uZ,3(τ )dτ " (17)
Define the parameter vector:

Λ = λ 1 λ 2 λ 3 T .
Equations ( 15), ( 16) and (17) allow to obtain the expressions of the estimate of the parameters λ1 , λ2 et λ3 :

D(t) Λ = N (t) (18)
with the matrices

D(t) ∈ R 3×3 et N (t) ∈ R 3×1 . If D(t) is invertible, one gets:   λ1 λ2 λ3   = D(t) -1 N (t) (19)
Remark 1: Note that at time t = 0, the matrices and vectors used to compute the estimations are null. The parameters are then undetermined. Therefore, the formula has to be calculated, not at the time t = 0, but at a later time, t = ǫ with ǫ > 0 and small. Thus the parameter vector is estimated as follows :

Λ = arbitrary for t ∈ [0, ǫ[ D(t) -1 N (t) for t ∈ [ǫ, ∞[

B. Implementation

In the previous paragraph, it was assumed that the perturbations were constant. Nevertheless, this is physically not realistic. It has been seen in Section II that there are two different perturbation sources. The first perturbation source, due to simplifications and modeling error, can be considered as a constant if the signals are integrating over a small sliding window. So far, the integrals have been computed over an interval [0, t]. However integrating over a short time interval [0, T f ] (T f represents the window length), we can consider that the perturbations are constant. Finally we use a change of variable to reduce the estimation interval

I T f 0 = [0, T f ] to [0, 1]. N (t)
and D(t) can be rewritten as:

N (t) = 1 0 g z,i (τ )Z(τ )dτ 1≤i≤3 and D(t) = 1 0 g uz,i (τ )U Z,j (τ )dτ 1≤i≤3,1≤j≤3 with g z,1 (τ ) = ˆ6 3! (1 -τ ) 3 -18 2 (1 -τ ) 2 τ + 9(1 -τ )τ 2 -τ 3 ˜T 3 f g z,2 (τ ) = ˆ-24 3! (1 -τ ) 3 τ + 36 2 (1 -τ ) 2 τ 2 -12(1 -τ )τ 3 + τ 4 ˜T 4 f g z,3 (τ ) = ˆ60 3! (1 -τ ) 3 τ 2 -60 2 (1 -τ ) 2 τ 3 + 15(1 -τ )τ 5 -τ 5 ˜T 5 f gu z,1 (τ ) = ˆ3 3! (1 -τ ) 3 τ 2 -1 2 (1 -τ ) 2 τ 3 ˜T 5 f gu z,2 (τ ) = ˆ-4 3! (1 -τ ) 3 τ 3 + 1 2 (1 -τ ) 2 τ 4 ˜T 6 f gu z,3 (τ ) = ˆ5 3! (1 -τ ) 3 τ 4 -1 2 (1 -τ ) 2 τ 5 ˜T 7 f
Using the trapezoidal method, the integrals are obtained from the output of a classical FIR filter. To reject the harmonic perturbations with frequencies close to rotational speed, the resulting numerator and denominator of each estimated parameters are filtered using a low-pass filter. The filter used is similar to the one presented in [START_REF] Trapero | An algebraic frequency estimator for a biased and noisy sinusoidal signal[END_REF] :

λi,f (t) = F (s)n(t) F (s)d(t) (20)
where F (s) =

IV. RESULTS

Simulation and experimental results are presented in this section. The numerical simulations are carried out in order to show the performances of the algorithm. Finally, we give the estimation results on the real process. The parameters are the rotor mass m = 6.7(kg) and the nominal air gap s = 5.10 -4 (m). The control is realized as a cascade, with a current controller in the inner loop and a position controller in the outer loop. The inner loop can be characterised by a control based on the electrical model of the bearing coils. The outer is based on the rigid body 6)-( 7). The chosen control design is a flatness based trajectory tracking control, as described in [START_REF] Lowis | Flatness-based trajectory tracking control of a rotationg shaft[END_REF]. The reference trajectory is an elliptic one where Y * = r Y cos(ωt) and Z * = r Z sin(ωt) with r Y = 30.10 -6 (m) and r Z = 70.10 -6 (m). The shaft is rotating with an angular velocity of about 3000 rpm, i.e. ω = 50 • 2π(rad/s). The estimator is implemented in a discretized way with a sampling frequency of 1.10 -4 (s). The choice of the window length is not straightforward. The window has to be sufficiently long to cancel the noise but not too long to keep the assumption on the constant perturbation true. The experience shows that a window T f = 0.02s gives good results. The initial value of λ k,0 is set arbitrarily. The filter F (s) is chosen with ζ = 0.707 and w n = 15rad/s. At the beginning of the experiment, the controller uses the initial and arbitrary λ k,0 from t = 0 to t = t e , t e being the convergence time of the algorithm. The estimator, connected in parallel to the system, estimates the parameters λk in real time on the interval (ǫ, t e ]. As soon as the estimations have converged to constant values, the initial values λ k,0 are replaced by the estimated ones.

A. Simulation Results

In the simulations, the real values of the parameters λ k are chosen as λ 1 = 1.8.10 -6 , λ 2 = 2.0.10 -6 , λ 1 = 2.2.10 -6 . The initial value of the estimation is λ k,0 = 4.0.10 -6 . To estimate the parameters, we choose sinusoidal perturbations with non zero mean value to be close to the real perturbations. Fig. 2 presents the parameters estimation. In this figure, it can be seen that the three parameters converge to the desired value in less than 0.1 second with λ1 = 1.786.10 -6 , λ2 = 2.072.10 -6 and λ3 = 2.199.10 -6 . It represents an error less that 3.5%. Fig. 3 depicts the trajectory tracking evolution of the z-axis only (the y-axis is similar). The estimated parameters are updated after 0.2 seconds. Before identification, the position errors are already small. However, the perturbations estimation is far from the real perturbation. It reflects the modeling errors. After parameters update, the position error is slightly decreased while the perturbation error is considerably decrease. The model is more correct, and the observed perturbation is close to the real one.

B. Experimental Results

In this Section, the experimental results are presented. The test-bench is shown Fig. 4. The computer hardware on the All currents are generated by 3 Dc brush amplifiers that serve as three independent control inputs (Fig. 5). When the experiment starts, one does not have precise knowledge of λ k . As in the simulations, the parameters are set arbitrarily equal to 4.0.10 -6 (N m 2 /A 2 ). The estimation algorithm starts at the beginning of the experiment. Fig. 6 shows that after 0.2 seconds, the parameters have converged. Indeed, one obtains three constant values for the parameters ( λ1 = 2.53.10 -6 , λ2 = 2.04.10 -6 and λ1 = 2.42.10 -6 ).

In order to check estimation accuracy of the parameters, we observe the trajectories and the perturbations before and after estimation. Fig. 7 and8 show the evolution of the y and z axes before identification while Fig. 9 and 10 show the evolution after identification. Similarly to the simulation results, the most obvious result is on the perturbations. An estimate of the latter is obtained using the perturbation observer given in [START_REF] Lowis | Flatness-based trajectory tracking control of a rotationg shaft[END_REF]. The amplitude of the perturbations is Fig. 8. Trajectory tracking evaluation on the z axe before identification at least three times smaller after identification than before and leads to a significant reduction of the modeling uncertainties. Concerning the trajectories tracking and their errors, the improvement is less sensitive. Looking at the errors, it can be seen that the positions errors are slightly lower when the parameters are well estimated and more regular. The behavior is not really different, but this is not surprising. Indeed, the designed control law compensate for the perturbations which are themselves estimated by an observer. Estimating the perturbation, the observer takes in consideration the modelling error. Finally the last plots (Fig. 11 and12) gives an idea of the shape of the hole before and after identification.

V. CONCLUSIONS AND FUTURE WORKS

In this article, it has been presented an estimator able to evaluate in a fast way and simultaneously several unknown parameters of a radial magnetic bearing. This system is strongly nonlinear and unstable. The performances of this algorithm show good experimental results. Indeed, the estimation is robust, the estimation error is low in spite of noise and perturbations, and also fast with convergence time of 0.2 seconds. Using the estimated parameters, the trajectories are well tracked. Experiments show that the maximum tracking error is around 5µm for a non-circular hole. As improvement, it will by interesting to simultaneously identify the perturbations amplitude. Using better sensors than the ones used could also leads to better tracking performance.

ω 2 n s 2 +2ζωn+ω 2 n.

 22 Experimental results are shown in the next Section.

2 )Fig. 2 .

 22 Fig. 2. Parameter estimation without noise

Fig. 3 .Fig. 4 .

 34 Fig. 3. Trajectory tracking evolution of z axe

  Fig. 5. Control Currents

Fig. 10 .

 10 Fig. 10. Trajectory tracking evaluation on the z axe after identification

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Position in the yz plan before identification

  Fig.9. Trajectory tracking evaluation on the y axe after identification

	z(t) (µ.m)	0 50										z * (t ) z (t )
		-50
			0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
								time (s)
		10
	ez(t )(µ . m)	-5 0 5										ez(t )
		-10	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
								time (s)
		10										p z( t )
	p z (t) (m)	-5 0 5
		-10
			0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
								time (s)
													40
													y * (t )
													y(t) (µ.m)	-20 0 20	y (t )
													-40	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
													time (s)
													10
													ey(t )(µ . m)	-5 0 5	ey(t )
													-10	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
													time (s)
													5
													p y( t )
													p y (t) (m)	-5 0
													-10	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
													time (s)

This work was supported by Nord-Pas de Calais Regional Council and FEDER through the Contrat de Projets Etat Region (CPER) 2007-2013.


Emmanuel Morin
email: emmanuel.morin@univ-nantes.fr

Emmanuel Prochasson
email: emmanuel.prochasson@univ-nantes.fr

Bilingual Lexicon Extraction from Comparable Corpora Enhanced with Parallel Corpora

In this article, we present a simple and effective approach for extracting bilingual lexicon from comparable corpora enhanced with parallel corpora. We make use of structural characteristics of the documents comprising the comparable corpus to extract parallel sentences with a high degree of quality. We then use state-of-the-art techniques to build a specialized bilingual lexicon from these sentences and evaluate the contribution of this lexicon when added to the comparable corpus-based alignment technique. Finally, the value of this approach is demonstrated by the improvement of translation accuracy for medical words.

Introduction

Bilingual lexicons are important resources of many applications of natural language processing such as cross-language information retrieval or machine translation. These lexicons are traditionally extracted from bilingual corpora.

In this area, the main work involves parallel corpora, i.e. a corpus that contains source texts and their translations. From sentence-to-sentence aligned corpora, symbolic [START_REF] Carl | An Intelligent Terminology Database as a Pre-processor for Statistical Machine Translation[END_REF], statistical [START_REF] Daille | Towards Automatic Extraction of Monolingual and Bilingual Terminology[END_REF], or hybrid techniques [START_REF] Gaussier | Modèles statistiques pour l'extraction de lexiques bilingues[END_REF] are used for word and expression alignments. However, despite good results in the compilation of bilingual lexicons, parallel corpora are rather scarce resources, especially for technical domains and for language pairs not involving English. For instance, current resources of parallel corpora are built from the proceedings of international institutions such as the European Union (11 languages) or the United Nations (6 languages), bilingual countries such as Canada (English and French languages), or bilingual regions such as Hong Kong (Chinese and English languages).

For these reasons, research in bilingual lexicon extraction is focused on another kind of bilingual corpora. These corpora, known as comparable corpora, are comprised of texts sharing common features such as domain, genre, register, sampling period, etc. without having a source text-target text relationship. Although the building of comparable corpora is easier than the building of parallel corpora, the results obtained thus far on comparable corpora are contrasted. For instance, good results are obtained from large corpora -several million words -for which the accuracy of the proposed translation is between 76% [START_REF] Fung | A Statistical View on Bilingual Lexicon Extraction: From Parallel Corpora to Non-parallel Corpora[END_REF] and 89% [START_REF] Rapp | Automatic Identification of Word Translations from Unrelated English and German Corpora[END_REF] for the first 20 candidates. [START_REF] Cao | Base Noun Phrase Translation Using Web Data and the EM Algorithm[END_REF] have achieved 91% accuracy for the top three candidates using the Web as a comparable corpus. But for technical domains, for which large corpora are not available, the results obtained, even though encouraging, are not completely satisfactory yet. For instance, [START_REF] Déjean | Une nouvelle approche à l'extraction de lexiques bilingues à partir de corpus comparables[END_REF] obtained a precision of 44% and 57% for the first 10 and 20 candidates in a 100,000-word medical corpus, and 35% and 42% in a multi-domain 8 million-word corpus. For French/English single words, [START_REF] Chiao | Looking for candidate translational equivalents in specialized, comparable corpora[END_REF] using a medical corpus of 1.2 million words, obtained a precision of about 50% and 60% for the top 10 and top 20 candidates. [START_REF] Morin | Bilingual Terminology Mining -Using Brain, not brawn comparable corpora[END_REF] obtained a precision of 51% and 60% for the top 10 and 20 candidates in a 1.5 million-word French-Japanese diabetes corpus.

The above work in bilingual lexicon extraction from comparable corpora relies on the assumption that words which have the same meaning in different languages tend to appear in the same lexical contexts [START_REF] Fung | A Statistical View on Bilingual Lexicon Extraction: From Parallel Corpora to Non-parallel Corpora[END_REF][START_REF] Rapp | Automatic Identification of Word Translations from Unrelated English and German Corpora[END_REF]. Based on this assumption, a standard approach consists of building context vectors for each word of the source and target languages. The candidate translations for a particular word are obtained by comparing the translated source context vector with all target context vectors. In this approach, the translation of the words of the source context vectors depends on the coverage of the bilingual dictionary vis-à-vis the corpus. This aspect can be a potential problem if too few corpus words are found in the bilingual dictionary [START_REF] Chiao | The effect of a general lexicon in corpus-based identification of French-English medical word translations[END_REF][START_REF] Déjean | Une nouvelle approche à l'extraction de lexiques bilingues à partir de corpus comparables[END_REF].

In this article, we want to show how this problem can be partially circumvented by combining a general bilingual dictionary with a specialized bilingual dictionary based on a parallel corpus extracted through mining of the comparable corpus. In the same way that recent works in Statistical Machine Translation (SMT) mines comparable corpora to discover parallel sentences [START_REF] Resnik | The Web as a parallel corpus[END_REF][START_REF] Yang | Automatic construction of English/Chinese parallel corpora[END_REF][START_REF] Dragos | Improving Machine Translation Performance by Exploiting Non-Parallel Corpora[END_REF]Abdul-Rauf and Schwenk, 2009, among others), this work contributes to the bridging of the gap between comparable and parallel corpora by offering a framework for bilingual lexicon extraction from comparable corpus with the help of parallel corpusbased pairs of terms.

The remainder of this article is organized as follows. In Section 2, we first present the method for bilingual lexicon extraction from comparable corpora enhanced with parallel corpora and the associated system architecture. We then quantify and analyse in Section 3 the performance improvement of our method on a medical comparable corpora when used to extract specialized bilingual lexicon. Finally, in Section 4, we discuss the present study and present our conclusions.

System Architecture

The overall architecture of the system for lexical alignment is shown in Figure 1 and comprises parallel corpus-and comparable corpus-based align-ments. Starting from a comparable corpus harvested from the web, we first propose to extract parallel sentences based on the structural characteristics of the documents harvested. These parallel sentences are then used to build a bilingual lexicon through a tool dedicated to bilingual lexicon extraction. Finally, this bilingual lexicon is used to perform the comparable corpus-based alignment. For a word to be translated, the output of the system is a ranked list of candidate translations.

Extracting Parallel Sentences from Comparable Corpora

Parallel sentence extraction from comparable corpora has been studied by a number of researchers [START_REF] Ma | Bits: A Method for Bilingual Text Search over the Web[END_REF][START_REF] Chen | Parallel Web Text Mining for Cross-Language Information Retrieval[END_REF][START_REF] Resnik | The Web as a parallel corpus[END_REF][START_REF] Yang | Automatic construction of English/Chinese parallel corpora[END_REF][START_REF] Fung | Mining Very-Non-Parallel Corpora: Parallel Sentence and Lexicon Extraction via Bootstrapping and EM[END_REF][START_REF] Dragos | Improving Machine Translation Performance by Exploiting Non-Parallel Corpora[END_REF][START_REF] Abdul | On the use of comparable corpora to improve SMT performance[END_REF], among others) and several systems have been developed such as BITS (Bilingual Internet Test Search) [START_REF] Ma | Bits: A Method for Bilingual Text Search over the Web[END_REF], PTMiner (Parallel Text Miner) [START_REF] Chen | Parallel Web Text Mining for Cross-Language Information Retrieval[END_REF], and STRAND (Structural Translation Recognition for Acquiring Natural Data) [START_REF] Resnik | The Web as a parallel corpus[END_REF]. Their work relies on the observation that a collection of texts in different languages composed independently and based on sharing common features such as content, domain, genre, register, sampling period, etc. contains probably some sentences with a source text-target text relationship. Based on this observation, dynamic programming [START_REF] Yang | Automatic construction of English/Chinese parallel corpora[END_REF], similarity measures such as Cosine [START_REF] Fung | Mining Very-Non-Parallel Corpora: Parallel Sentence and Lexicon Extraction via Bootstrapping and EM[END_REF] or word and translation error ratios [START_REF] Abdul | On the use of comparable corpora to improve SMT performance[END_REF], or maximum entropy classifier [START_REF] Dragos | Improving Machine Translation Performance by Exploiting Non-Parallel Corpora[END_REF] are used for discovering parallel sentences.

Although our purpose is similar to these works, the amount of data required by these techniques makes them ineffective when applied to specialized comparable corpora used to discover parallel sentences. In addition, the focus of this paper is not to propose a new technique for this task but to study how parallel sentences extracted from a comparable corpus can improve the quality of the candidate translations. For theses reasons, we propose to make use of structural characteristics of the documents comprising the comparable corpus to extract auto- In fact, specialized comparable corpora are generally constructed via the consultation of specialized Web portals. For instance, [START_REF] Chiao | Looking for candidate translational equivalents in specialized, comparable corpora[END_REF] use CISMeF1 for building the French part of their comparable corpora and CliniWeb2 for the English part, and [START_REF] Déjean | Une nouvelle approche à l'extraction de lexiques bilingues à partir de corpus comparables[END_REF] use documents extracted from MEDLINE3 to build a German/English comparable corpus. Consequently, the documents collected through these portals are often scientific papers. Moreover, when the language of these papers is not the English, the paper usually comprises an abstract, keywords and title in the native language and their translations in the English language. These characteristics of scientific paper is useful for the efficient extraction of parallel sentences or word translations from the documents forming a specialized comparable corpus for which one part will inevitably be in English.

In this study, the documents comprising the French/English specialized comparable corpus were taken from the medical domain within the subdomain of 'breast cancer'. These documents have been automatically selected from the Elsevier website4 among the articles published between 2001 and 2008 for which the title or the keywords of the articles contain the multi-word term 'cancer du sein' in French and 'breast cancer' in English. We thus collected 130 documents in French and 118 in English and about 530,000 words for each language. Since the 130 French documents previously collected are scientific papers, each document contains a French abstract which is accompanied by its English translation. We exploit this structural characteristic of the French documents in order to build a small specialized parallel corpus directly correlated to the subdomain of 'breast cancer' involved in the comparable corpus.

Parallel Corpus-Based Alignment

We use the Uplug5 collection of tools for alignment [START_REF] Tiedemann | Recycling Translations -Extraction of Lexical Data from Parallel Corpora and their Application in Natural Language Processing[END_REF] to extract translations from our specialized parallel corpus. The output of such a tool is a list of aligned parts of sentences, that has to be post-process and filtered in our case. We clean the alignment with a simple yet efficient method in order to obtain only word translations. We associate every source word from a source sequence with every target word from the target sequence. As an example, uplug efficiently aligns the English word breast cancer with the French word cancer du sein (the data are described in Section 3.1). We obtain the following lexical alignment:

• cancer (fr) → (en) breast, cancer

• du (fr) → (en) breast, cancer

• sein (fr) → (en) breast, cancer
With more occurrences of the French word cancer, we are able to align it with the English words {breast, cancer, cancer, cancer, the, of, breast, can-cer}. We can then filter such a list by counting the translation candidates. In the previous example, we obtain: cancer (fr) → breast/2, the /1, of/1, cancer/4. The English word cancer is here the best match for the French word cancer. In many cases, only one alignment is obtained. For example, there is only one occurrence of the French word chromosome, aligned with the English word chromosome.

In order to filter translation candidates, we keep 1:1 candidates if their frequencies are comparable in the original corpus. We keep the most frequent translation candidates (in the previous example, cancer) if their frequencies in the corpus are also comparable. This in-corpus frequency constraint is useful for discarding candidates that appear in many alignments (such as functional words). The criterion for frequency acceptability is:

min(f 1 , f 2 )/max(f 1 , f 2 ) > 2/3
with f 1 and f 2 the frequency of words to be aligned in the parallel corpus.

By this way, we build a French/English specialized bilingual lexicon from the parallel corpus. This lexicon, called breast cancer dictionary (BC dictionary) in the remainder of this article, is composed of 549 French/English single words.

Comparable Corpus-Based Alignment

The comparable corpus-based alignment relies on the simple observation that a word and its translation tend to appear in the same lexical contexts. Based on this observation, the alignment method, known as the standard approach, builds context vectors in the source and the target languages where each vector element represents a word which occurs within the window of the word to be translated (for instance a seven-word window approximates syntactic dependencies). In order to emphasize significant words in the context vector and to reduce word-frequency effects, the context vectors are normalized according to association measures. Then, the translation is obtained by comparing the source context vector to each translation candidate vector after having translated each element of the source vector with a general dictionary.

The implementation of this approach can be carried out by applying the four following steps [START_REF] Fung | A Statistical View on Bilingual Lexicon Extraction: From Parallel Corpora to Non-parallel Corpora[END_REF][START_REF] Rapp | Automatic Identification of Word Translations from Unrelated English and German Corpora[END_REF]:

1. We collect all the lexical units in the context of each lexical unit i and count their occurrence frequency in a window of n words around i.

For each lexical unit i of the source and the target languages, we obtain a context vector v i which gathers the set of co-occurrence units j associated with the number of times that j and i occur together occ(i, j). In order to identify specific words in the lexical context and to reduce word-frequency effects, we normalize context vectors using an association score such as Mutual Information (MI) or Log-likelihood, as shown in equations 1 and 2 and in Table 1 (where

N = a + b + c + d).
2. Using a bilingual dictionary, we translate the lexical units of the source context vector. If the bilingual dictionary provides several translations for a lexical unit, we consider all of them but weight the different translations according to their frequency in the target language.

3. For a lexical unit to be translated, we compute the similarity between the translated context vector and all target vectors through vector distance measures such as Cosine or Weighted Jaccard (WJ) (see equations 3 and 4 where assoc i j stands for "association score").

4. The candidate translations of a lexical unit are the target lexical units ranked following the similarity score.

j ¬j i a = occ(i, j) b = occ(i, ¬j) ¬i c = occ(¬i, j) d = occ(¬i, ¬j) Table 1: Contingency table M I(i, j) = log a (a + b)(a + c) (1) λ(i, j) = a log(a) + b log(b) + c log(c) +d log(d) + (N ) log(N ) -(a + b) log(a + b) -(a + c) log(a + c) -(b + d) log(b + d) -(c + d) log(c + d) (2) Cosine v k v l = t assoc l t assoc k t t assoc l t 2 t assoc k t 2
(3)

W J v k v l = t min(assoc l t , assoc k t ) t max(assoc l t , assoc k t ) (4)
This approach is sensitive to the choice of parameters such as the size of the context, the choice of the association and similarity measures. The most complete study about the influence of these parameters on the quality of bilingual alignment has been carried out by [START_REF] Laroche | Revisiting Context-based Projection Methods for Term-Translation Spotting in Comparable Corpora[END_REF].

Experiments and Results

In the previous section, we have introduced our comparable corpus and described the method dedicated to bilingual lexicon extraction from comparable corpora enhanced with parallel corpora. In this section, we then quantify and analyse the performance improvement of our method on a medical comparable corpus when used to extract specialized bilingual lexicon.

Experimental Test bed

The documents comprising the French/English specialized comparable corpus have been normalised through the following linguistic pre-processing steps: tokenisation, part-of-speech tagging, and lemmatisation. Next, the function words were removed and the words occurring less than twice in the French and the English parts were discarded. Finally, the comparable corpus comprised about 7,400 distinct words in French and 8,200 in English.

In this study, we used four types of bilingual dictionary: i) the Wiktionary6 free-content multilingual dictionary, ii) the ELRA-M00337 professional French/English bilingual dictionary, iii) the MeSH8 metha-thesaurus, and iv) the BC dictionary (see Section 2.2). Table 2 shows the main features of the dictionaries, namely: the number of distinct French single words in the dictionary (# SWs dico.), the number of distinct French single words in the dictionary after projection on the French part of the comparable corpus (# SWs corpus), and the number of translations per entry in the dictionary (# TPE). For instance, 42% of the French context vectors could be translated with the Wiktionary (3,099/7,400). In bilingual terminology extraction from specialized comparable corpora, the terminology reference list required to evaluate the performance of the alignment programs are often composed of 100 single-word terms (SWTs) (180 SWTs in [START_REF] Déjean | Une nouvelle approche à l'extraction de lexiques bilingues à partir de corpus comparables[END_REF], 95 SWTs in [START_REF] Chiao | Looking for candidate translational equivalents in specialized, comparable corpora[END_REF], and 100 SWTs in [START_REF] Daille | French-English Terminology Extraction from Comparable Corpora[END_REF]). To build our reference list, we selected 400 French/English SWTs from the UMLS 9 meta-thesaurus and the Grand dictionnaire terminologique 10 . We kept only the French/English pair of SWTs which occur more than five times in each part of the comparable corpus. As a result of filtering, 122 French/English SWTs were extracted.

Experimental Results

In order to evaluate the influence of the parallel corpus-based bilingual lexicon induced from the comparable corpus on the quality of comparable corpus based-bilingual terminology extraction, four experiments were carried out. For each experiment, we change the bilingual dictionary required for the translation phase of the standard approach (see Section 2.3):

1. The first experiment uses only the Wiktionary.

Since the coverage of the Wiktionary from the comparable corpus is small (see Table 2), the results obtained with this dictionary yield a lower boundary.

2. The second experiment uses the Wiktionary added to the BC dictionary. This experiment attempts to verify the hypothesis of this study.

3. The third experiment uses the Wiktionary added to the MeSH thesaurus. This experiment attempts to determine whether a specialised dictionary (in this case the MeSH) would be more suitable than a specialized bilingual dictionary (in this case the BC dictionary) directly extracted from the corpus.

4. The last experiment uses only the ELRA dictionary. Since the coverage of the ELRA dictionary from the comparable corpus is the best (see Table 2), the results obtained with this one yield a higher boundary.

Table 3 shows the coverage of the four bilingual lexical resources involved in the previous experiments in the comparable corpus. The first column indicates the number of single words belonging to a dictionary found in the comparable corpus (# SWs corpus). The other column indicates the coverage of each dictionary in the ELRA dictionary (Coverage ELRA). Here, 98.9% of the single words belonging to the Wiktionary are included 10 http://www.granddictionnaire.com/ in the ELRA dictionary whereas less than 95% of the single words belonging to the Wiktionary+BC and Wiktionary+MeSH dictionaries are included in the ELRA dictionary. Moreover, the MeSH and BC dictionaries are two rather distinct specialized resources since they have only 117 single words in common. In the experiments reported here, the size of the context window n was set to 3 (i.e. a seven-word window), the association measure was the Mutual Information and the distance measure the Cosine (see Section 2.3). Other combinations of parameters were assessed but the previous parameters turned out to give the best performance.

Figure 2 summarises the results obtained for the four experiments for the terms belonging to the reference list according to the French to English direction. As one could expect, the precision of the result obtained with the ELRA dictionary is the best and the precision obtained with the Wiktionary is the lowest. For instance, the ELRA dictionary improves the precision of the Wiktionary by about 14 points for the Top 10 and 9 points for the top 20. These results confirm that the coverage of the dictionary is an important factor in the quality of the results obtained. Now, when you add the BC dictionary to the Wiktionary, the results obtained are also much better than those obtained with the Wiktionary alone and very similar to those obtained with the ELRA dictionary alone (without taking into account the top 5). This result suggests that a standard general language dictionary enriched with a small specialized dictionary can replace a large general language dictionary.

Furthermore, this combination is more interesting than the combination of the MeSH dictionary with the Wiktionary. Since the BC dictionary is induced from the corpus, this dictionary is directly correlated to the theme of breast cancer involved in the corpus. Consequently the BC dictionary is more suitable than the MeSH dictionary i) even if the MeSH dictionary specializes in the medical domain and ii) even if more words in the comparable corpus are found in the MeSH dictionary than in the BC dictionary.

This last observation should make us relativize the claim: the greater the number of context vector elements that are translated, the more discriminating the context vector will be for selecting translations in the target language. We must also take into account the specificity of the context vector elements in accordance with the thematic of the documents making up the corpus studied in order to improve bilingual lexicon extraction from specialized comparable corpora.

Conclusion and Discussion

In this article, we have shown how the quality of bilingual lexicon extraction from comparable corpora could be improved with a small specialized bilingual lexicon induced through parallel sentences included in the comparable corpus. We have evaluated the performance improvement of our method on a French/English comparable corpus within the sub-domain of breast cancer in the medical domain. Our experimental results show that this simple bilingual lexicon, when combined with a general dictionary, helps improve the accuracy of single word alignments by about 14 points for the Top 10 and 9 points for the top 20. Even though we focus here on one structural characteristic (i.e. the abstracts) of the documents comprising the comparable corpus to discover parallel sentences and induced bilingual lexicon, the method could be easily applied to other comparable corpora for which a bilingual dictionary can be extracted by using other characteristics such as the presence of parallel segments or paraphrases in the documents making up the comparable corpus.

Figure 1 :

 1 Figure 1: Overview of the system for lexical alignment

Figure 2 :

 2 Figure 2: Precision of translations found according to the rank

Table 2 :

 2 Main features of the French/English dictionaries

	Name	# SWs # SWs # TPE
		dict. corpus
	Wiktionary 20,317	3,099	1.8
	ELRA	50,330	4,567	2.8
	MeSH	18,972	833	1.6
	BC	549	549	1.0

Table 3 :

 3 Coverage of the bilingual lexical resources in the comparable corpus

	Name	# SWs Coverage
		corpus	ELRA
	Wiktionary	3,099	98.8%
	Wiktionary + BC	3,326	94.8%
	Wiktionary + MeSH	3,465	94.9%
	ELRA	4,567	100%

Proceedings of the 4th Workshop on Building and Using Comparable Corpora, pages 27-34, 49th Annual Meeting of the Association for Computational Linguistics, Portland, Oregon, 24 June 2011. c 2011 Association for Computational Linguistics

http://www.chu-rouen.fr/cismef/

http://www.ohsu.edu/cliniweb/

http://www.ncbi.nlm.nih.gov/PubMed

http://www.elsevier.com

http://stp.ling.uu.se/cgi-bin/joerg/ Uplug

http://www.wiktionary.org/

http://www.elra.info/

http://www.ncbi.nlm.nih.gov/mesh

http://www.nlm.nih.gov/research/umls

Acknowledgments

The research leading to these results has received funding from the European Communitys Seventh Framework Programme (FP7/2007-2013) under Grant Agreement no 248005.