Jean-François Patri

Julien Diard

Pascal Perrier

Optimal speech motor control and token-to-token variability: A Bayesian modeling approach

Keywords: Speech motor control, Speech sequence motor planning, Bayesian modeling, Optimal motor control

The remarkable capacity of the speech motor system to adapt to various speech conditions is due to an excess of degrees of freedom, which enables producing similar acoustical properties with different sets of control strategies. To explain how the Central Nervous System selects one of the possible strategies, a common approach, in line with optimal motor control theories, is to model speech motor planning as the solution of an optimality problem based on cost functions. Despite the success of this approach, one of its drawbacks is the intrinsic contradiction between the concept of optimality and the observed experimental intra-speaker token-to-token variability. The present paper proposes an alternative approach by formulating feedforward optimal control in a probabilistic Bayesian modeling framework. This is illustrated by controlling a biomechanical model of the vocal tract for speech production and by comparing it with an existing optimal

Introduction

Motor control aims at finding patterns of activation that agents should generate in their articulatory chain in order to achieve desired motor goals. This is in essence an ill-posed problem, since degrees of freedom of articulatory chains often largely exceed the degrees of freedom of the task. Therefore there is a multiplicity of possible solutions for achieving the desired motor goal. Optimal motor control theories aim at resolving this well-known redundancy problem [START_REF] Jordan | Computational motor control[END_REF][START_REF] Uno | Formation control of optimal trajectory in human multijoint arm movement: Minimum torque-change model[END_REF] by considering a cost function that attributes to each possible solution a certain performance value. The redundancy is resolved by this mean, if there is a unique solution that optimizes the value associated with this criterion.

A crucial consequence of this approach is that the resulting behavior of the controlled system is stereotyped. This means that for a given task and in a specified condition, the optimal solution is always the same and no trial-to-trial variability can be obtained. While this may be desirable in engineering applications, it is a major drawback for models that try to reproduce behavior of biological agents. In that sense, most optimal motor control models have been successful in accounting for average patterns of behavior at the expense of not being able to model trial-to-trial variability (see [START_REF] Todorov | Optimality principles in sensorimotor control[END_REF] for a review). Even though this issue has been addressed by stochastic optimal feedback control [START_REF] Todorov | Optimal feedback control as a theory of motor coordination[END_REF] by explicitly considering feedback in the planning process (leading to a closed loop control), this approach only concerns movements that can use on-line feedback information. In the case of very fast movements though, ongoing control is unlikely to rely on feedback, due to delays in afferent signals (for example 30-100 ms for visual motion in Schmolesky et al (1998); 99-143 ms for auditory feedback in speech processing in [START_REF] Tourville | Neural mechanisms underlying auditory feedback control of speech[END_REF]), and control is rather assumed to be performed through an open loop planning (see for example [START_REF] Kawato | Internal models for motor control and trajectory planning[END_REF]).

Nevertheless, two approaches are usually considered for recovering variability under a feedforward optimal control model. The first approach assumes that optimal planning is stereotyped, but that variability arises from noise in the pathway of the control signal or in the dynamics of the articulatory chain. This is the approach followed by stochastic optimal control theory. The second approach assumes that planning is driven by optimality, but that its realization does not systematically lead to the unique optimal situation. It should be noted that even if the control process is certainly subject to stochastic dynamics, the first approach alone is not satisfactory for explaining situations where variability results from different specific patterns of behavior (see for example [START_REF] Shim | Prehension synergies: trial-to-trial variability and hierarchical organization of stable performance[END_REF] for prehension tasks). Rather, the second approach better accounts for systematic deviations from a single optimal solution. For instance the role of motor memory on convergence to locally optimal solutions rather than a global optimum has been suggested as a crucial aspect of motor control [START_REF] Ganesh | Motor memory and local minimization of error and effort, not global optimization, determine motor behavior[END_REF]. These questions are of particular interest in speech motor control. Indeed, it is unlikely that speech control relies primarily on feedback signals, due to the speed of tongue movement [START_REF] Perkell | Speech motor control: Acoustic goals, saturation effects, auditory feedback and internal models[END_REF]). Yet, trial-totrial variability is observed in phoneme production at the acoustic, articulatory and muscle activation levels [START_REF] Perkell | Variability in production of the vowels /i/ and /a[END_REF]. This variability is underpinned by the presence of redundancy at the three levels described above: 1) a particular phonemic goal does not correspond to a unique point in the acoustic domain since different acoustic signals are perceived as a unique phonemic category, 2) a particular acoustic signal can be produced by different vocal tract configurations, and 3) a particular vocal tract configuration can be attained by different patterns of muscle activation.

Attempts at modeling feedforward speech motor control based on optimal control theory have been able to reproduce a number of experimental speech patterns [START_REF] Ma | Anticipatory coarticulation in vowel-consonant-vowel sequences: A crosslinguistic study of french and mandarin speakers[END_REF][START_REF] Perrier | Modeling the production of VCV sequences via the inversion of a biomechanical model of the tongue[END_REF][START_REF] Guenther | A theoretical investigation of reference frames for the planning of speech movements[END_REF][START_REF] Guenther | Speech sound acquisition, coarticulation, and rate effects in a neural network model of speech production[END_REF]. However, while these results are consistent with average values among and across subjects, they fail at accounting for individual trial-to-trial variability from a theoretical point of view.

The present work aims at addressing this issue. By formulating optimal control in a Bayesian modeling framework we suggest that both variability and selection of motor control variables in speech production can be obtained in a principled way, from uncertainty at the representational level and without resorting solely to stochastic noise in the dynamics. We illustrate this approach by presenting a Bayesian formulation of an optimal control model for speech motor planning [START_REF] Perrier | Modeling the production of VCV sequences via the inversion of a biomechanical model of the tongue[END_REF].

The remainder of this paper is divided into four sections. Section 2 describes the essential ingredients of GEPPETO, the optimal control model that we aim to reformulate. From these ingredients the Bayesian model is then introduced in Section 3, in two steps. The first step consists in a Bayesian model inferring motor control variables for the production of a single phoneme. The second step consists in the complete Bayesian formulation of GEPPETO, that is, a Bayesian model planning optimal motor control variables for the production of sequences of phonemes. In Section 4 we compare and discuss the main results of the Bayesian formulation with respect to its optimal control version.

Main ingredients of GEPPETO

This section summarizes the key components of GEP-PETO, the optimal control model for speech motor planning that we aim at reformulating in the Bayesian framework. The following description focuses on the main hypotheses and we refer the reader to [START_REF] Perrier | Modeling the production of VCV sequences via the inversion of a biomechanical model of the tongue[END_REF] for a more detailed description of the model. We structure them as follows.

H 1 : GEPPETO computes the motor control variables of a biomechanical model of the tongue [START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the equilibrium point hypothesis[END_REF][START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF] in order to produce a desired sequence of speech gestures. GEPPETO defines a speech sequence as a succession of fundamental phonological units, corresponding to phonemes of the considered language. As the model only includes an account of the tongue, only phonemes that do not require lip rounding are considered. The set of considered phonemes is therefore {/i/, /e/, /E/, /a/, /oe/, /O/, /k/}. The first hypothesis of GEPPETO is therefore:

(H 1 ) The motor control of a speech sequence is organized on the basis of the specification of motor goals that are related to phonemes.

H 2 : GEPPETO supposes that phonemes are characterized and controlled in the acoustic domain and that the acoustic signal is characterized by the value of the first 3 peaks of the spectral envelope (these peaks are called "formants"). Therefore:

(H 2a ) The acoustic signal is represented by a point in a 3 dimensional space. (H 2b ) Phonemic goals are represented by specific simply connected regions of the 3-dimensional acoustic space. These regions are accounted for by ellipsoids defined through dispersion regions measured from phoneme production experiments [START_REF] Mt | Modèles d'intégration audiovisuelle de signaux linguistiques : de la perception humaine a la reconnaissance automatique des voyelles[END_REF][START_REF] Ménard | Production et perception des voyelles au cours de la croissance du conduit vocal : variabilité, invariance et normalisation[END_REF]Calliope, 1984). These regions are represented in Figure 1 by their projections on the (F 2 , F 1 ) and (F 2 , F 3 ) planes.

H 3 : The biomechanical model on which GEPPETO is based consists of a finite element structure representing the projection of the tongue on the mid-sagittal plane. Six principal muscles are considered as actuators for shaping the tongue. Figure 2 represents the tongue configuration at rest as well as the fibers of one of the considered muscles (the posterior genioglossus). The resulting tongue shape corresponds to the mechanical equilibrium of the forces generated by each muscle; the activation of each muscle is specified through a λ parameter, which specifies the muscle length above which active muscle force is generated 1 , in agreement with the Equilibrium Point Hypothesis [START_REF] Feldman | Once more on the equilibriumpoint hypothesis (λ model) for motor control[END_REF]. Hence:

(H 3 ) Achieving a particular configuration of the tongue consists in specifying a point in the 6-dimensional control space (λ 1 , . . . , λ 6 ).

1 Muscle force F generated by the biomechanical model is specified as

F = ρ[exp(cA) -1], ( 1
)
where c is a form parameter accounting for the gain of the feedback from the muscle to the motoneurons pool and ρ a magnitude parameter directly related to force-generating capability. A is the muscle activation corresponding to

A = l -λ + µ l, ( 2
)
where l is the actual muscle length, l the muscle shortening or lengthening velocity and µ a damping coefficient due to proprioceptive feedback [START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the equilibrium point hypothesis[END_REF].

H 4 : For every tongue configuration the resulting acoustic signal is generated from the computation of the vocal tract volume (via a model that links mid-sagittal views and cross-sectional areas from the glottis to the lips [START_REF] Perrier | Vocal tract area function estimation from midsagittal dimensions with ct scans and a vocal tract castmodeling the transition with two sets of coefficients[END_REF]). For every point in the 6dimensional control space there is therefore a unique associated point in the 3-dimensional acoustic space. Furthermore, GEPPETO assumes that:

(H 4 ) The knowledge of the mapping from the control variables to the acoustic domain is stored in an internal model in the Central Nervous System (CNS).

It is assumed that this model results from a learning process that generalizes the relation between motor control variables and formants from a limited number of examples. This model is considered to be "static" as it associates motor control variables and outputs at targets. It is implemented through a Radial Basis Function (RBF) network [START_REF] Poggio | A theory of networks for approximation and learning[END_REF]. This neural network is learned through classical supervised learning.

H 5 : GEPPETO assumes that articulatory trajectories between two successive targets emerge from the interactions between the motor control variables at targets, the specified duration of the transition between targets and the biomechanical properties of the tongue. It is important to note that GEPPETO does not assume any kind of specification of desired trajectories or any optimization of a cost at the level of the trajectories (such as maximum velocity, jerk, total amount of force). The specification of the appropriate control variables for the generation of a sequence does not involve inverse dynamics. Once the motor control variables at the targets are specified, the time variations of these variables are assumed to proceed from the value at target n to the values at target n + 1 at a constant rate of shift. This is consistent with the suggestion made by [START_REF] Laboissière | The control of multi-muscle systems: human jaw and hyoid movements[END_REF].

H 6 : The aim of GEPPETO is to specify in the 6dimensional control space a discrete sequence that generates a sequence of acoustic goals at the targets that are inside the ellipsoids characterizing the motor goals of the different phonemes. This is an ill-posed problem as there is an infinity of possible trajectories reaching the desired phonemic targets. To resolve this redundancy, GEPPETO assumes that the controller selects the trajectory that is optimal in the displacement of the corresponding λ variables (i.e. motor control space) .

To this end, GEPPETO defines a cost function measuring the distance between control variables over the whole set of targets within the sequence. For a threephoneme sequence, which will be our focus for the re-mainder of the text, this corresponds to the perimeter of the triangle defined by the three control points 2 . A subtlety about the interpretation of this cost function should be mentioned here. The cost function introduces a term relating the first and last phonemes in the sequence. This could be interpreted as sequence planning where only neighbors phonemes would influence each other, and where there would be a return, from the last phoneme, to the first one. However this is not what it is intended here. The term relating the first to the last phoneme is introduced in order to model forward and backward planned coarticulation influences. In this sense, in a sequence of 3 phonemes, this cost function introduces dependencies between items independently of their relative order in the sequence, allowing influence of every phoneme on every other one.

(H 6a ) For a three-phoneme sequence, the planning problem consists in finding a set of three points in the 6dimensional control space that minimizes the perimeter of the triangle they define, under the perceptual constraint that the corresponding spectral properties of the signal are within the ellipsoids regions assigned to each phoneme in the sequence. (H 6b ) The optimization process is performed by a gradient descent algorithm where the perceptual constraint is specified as an additional cost that vanishes whenever the corresponding acoustic signal falls within the correct ellipsoid region and goes to infinity (in practice, a large number) otherwise.

Bayesian formulation of GEPPETO

This section presents the Bayesian framework within which the reformulation of the optimal control model described in the previous section is derived. The approach is based on the Bayesian Programming methodology [START_REF] Lebeltel | Bayesian robot programming[END_REF]Bessière et al, 2013Bessière et al, , 2008) ) that proposes a precise structure for the definition of a Bayesian model. In order to illustrate the framework and to derive the general model stepwise, we first expose a sub-model, aiming at generating the control variables 2 For simplicity, the main text presents the case of sequences of 3-phonemes, without loss of generality. For a general n-phoneme sequence the proposed cost function would correspond to the perimeter of the corresponding (n -1)-simplex defined by the n control variables in the 6dimensional control space. For the present 3-phoneme case the 2-simplex corresponds to the triangle introduced in the text. Rigorously, influence of every phoneme of the sequence on every other one would be rather modeled by a cost function involving distances between every pair of phonemes. In order to avoid the corresponding quadratic combinatorial growth of the number of terms in the cost function, its definition has been simplified into the one presented here. for the production of a single phoneme, before focusing on the generation of sequences of phonemes.

3.1 Bayesian model for the production of a single phoneme

Description

This section presents the description of the Bayesian model, obtained by translating into probabilistic terms the hypotheses and knowledge introduced above.

Variables The variables of the Bayesian model are simply extracted from the key ingredients of GEPPETO described in Section 2.

Φ is the variable representing all the categories of phonemes.

It is a discrete variable composed by all the phonemes specified in hypothesis H 1 of Section 2. An additional "no-phoneme" category (denoted by /00/) is further assumed in order to take into account all acoustic configurations that do not fall within any of the above phonemic categories. The values taken by this variable are labeled by: Φ = { /i/, /e/, /E/, /a/, /oe/, /O/, /k/, /00/ } .

S represents the spectral characteristics of the acoustic signal. Hypothesis H 2a specifies this signal as a point in a 3-dimensional formant space. Therefore, S corresponds to a continuous vector variable, S = (F 1 , F 2 , F 3 ). The domain for this variable is the same as in GEPPETO and corresponds to the acoustic values attained by the simulations of the biomechanical model of the tongue. M represents the motor control variables controlling the articulatory configurations of the tongue. According to H 3 , these control variables correspond to the six λ parameters specifying the activation threshold length of each muscle. M is therefore a continuous 6-dimensional vector variable defined by M = (λ 1 , . . . , λ 6 ). The domain of M is specified as in GEPPETO and corresponds to the values of each λ for which the bio-mechanical model attains its equilibrium configurations, constrained by the vocal tract boundaries.

Decomposition We now define the structure of the Bayesian model, by specifying the joint probability distribution over the three above variables. Following the chain rule, an exact decomposition of the joint probability distribution P (M S Φ) is given by:

P (M S Φ) = P (M ) P (S | M ) P (Φ | S M ).
(3)

In order to avoid confusions, we draw attention to the notation that is employed here. The domain of the joint probability distribution constructed here is composed of discrete and continuous variables. Usually, one writes P for probability distributions over discrete variables and p for probability densities over continuous variables. For simplicity, we chose not to make this distinction here. Similarly, all summations and integrals are denoted by the sign , even when rigorously it is the sign that should be used for continous variables. Now, due to hypothesis H 2 described in Section 2, the last term of the decomposition of Equation ( 3) can be simplified. Indeed, according to this hypothesis phonemes are assumed to be fully characterized by their characteristics in acoustic space. Therefore, it can be assumed that Φ is independent of M conditioned on the knowledge of S. Under this assumption the joint probability distribution becomes

P (M S Φ) = P (M ) P (S | M ) P (Φ | S).
(4)

Figure 3 illustrates the Bayesian network representing this decomposition.

Parametric forms

P (M ) is the prior probability distribution over motor control variables M . Since no prior knowledge is assumed about this variable, P (M ) is defined as a uniform probability distribution over domain D M :

P (M ) = 1 |D M | if M ∈ D M 0 otherwise. (5)
P (S | M ) represents the knowledge about the spectral characteristics of the acoustic signal produced, given the equilibrium configuration of the tongue attained for the motor variables M . This knowledge is described in H 4 and corresponds to the internal model implemented by the RBF network of GEP-PETO. Denoting by S * (M ) the spectral properties of the acoustic signal associated to M by this RBF network, the corresponding probability distribution is assumed to be deterministic and is given by: Gaussian distribution is controlled by a parameter, denoted by κ S , that multiplies its variance. This allows to control the precision of the categorization task performed by the probability distribution P (Φ | S). This sub-model is not further described here and we refer the reader to the supplementary material for a more detailed description. For illustration, Figure 4 presents an example of probability distribution P (S | [Φ = k]) as well as the resulting likelihood function

P (S | M ) = δ S * (M ) (S) (6
P ([Φ = k] | S) = f (S).
Parameter identification The last step in order to completely specify the joint probability distribution given by Equation ( 4) is the identification of parameter values characterizing its probability distributions. κ S is the only parameter that remains unspecified. Its specification reflects the integration of particular assumptions in the model and its value will be given for each result. The effect of parameter κ S on the likelihood function P (Φ | S) is illustrated in Figure 5.

Inference of control variables M for the production of a given phoneme Φ

Having specified the joint probability distribution P (M S Φ), we now formulate the question to be solved by the Bayesian model. As the problem is to infer motor control variables M producing a desired phoneme Φ, the approach consists in computing the probability distribution over M , conditioned on the specified value of Φ.

The corresponding probability distribution, P (M | Φ), is obtained by standard Bayesian inference as

P (M | Φ) = P (M Φ) P (Φ) = S P (M S Φ) P (Φ) ∝ S P (M S Φ), (7)
where the proportionality symbol "∝" on the last line accounts for the term P (Φ), which does not depend on M for a given Φ value. Now, using the decomposition of Equation ( 4) we have:

P (M | Φ) ∝ S P (M S Φ) ∝ S P (M ) P (S | M ) P (Φ | S) ∝ S P (S | M )P (Φ | S) ∝ P (Φ | S * (M )), (8)
where the third line followed because P (M ) is assumed to be a uniform distribution and the last line is derived by recalling that P (S | M ) is zero unless S = S * (M ) (see Equation ( 6)).

Implementation of the model

The aim of the model is to generate motor control variables M performing a desired phoneme Φ. The probability distribution P (M | Φ) characterizes every control variable M with its probability for achieving the desired phoneme Φ. The best choice would be to select control variables that maximize this probability. However this would eliminate any possible variability and lead to the stereotyped situation encountered in GEPPETO. As the aim is to preserve variability, we adopted a decision policy based on a random sampling of the control variables space from P (M | Φ). Accuracy of the obtained solutions is ensured in average, since with this sampling the most probable control variables correspond to the ones having high probability of achieving the desired phoneme. This sampling is implemented by a standard Markov Chain Monte Carlo algorithm (MCMC), which performs a random walk that has the desired probability distribution as its equilibrium distribution. We draw attention on the interpretation of this particular implementation of the Bayesian model. We are not assuming that the biological system is indeed performing MCMC sampling. In terms of a biological implementation of this process, one would imagine that the brain stores information about P (M | Φ) in some ways and would use it to optimize the mapping from phoneme to motor space.

Results

Figure 6 shows histograms of control variables samples M , obtained from P (M | Φ) as described in the previous section. As they are 6-dimensional probability distributions, they are represented by their six marginal distributions. It can be noted that each phoneme corresponds to a specific set of distributions of control variables λ. Some of these control variables appear to be constrained within small ranges of values, for instance Small values of the λ control variables correspond to high levels of muscle activation and vice versa.

We can thus see that the Bayesian model correctly extracts the antagonist interaction of these two muscles in the front/high and back/low movement direction. This antagonism has been found in electromyographic measures of muscle activity during speech production [START_REF] Honda | Organization of tongue articulation for vowels[END_REF] Figure 2). This direction of movement is thus coherent with the variation of position of the tongue for the production of these four phonemes, qualitatively confirming the good adequacy of the Bayesian model with experimental results.

In order to assess the performance of the Bayesian model it is necessary to evaluate its capacity to effectively generate spectral properties that distribute around the correct areas defined in the acoustic space for each phoneme. Figure 7 represents the histograms of the first three formants of the acoustic signals corresponding to the samples M of Figure 6. These values were obtained through the RBF network that models the mapping between M and S, as described in Section 2. It can be seen that the obtained formants correctly distribute inside the goal regions.

Bayesian model for planning a sequence of phonemes

The previous section described a Bayesian model inferring motor control variables for the production of a unique phoneme. We now turn to planning sequences of phonemes under a "minimum effort" assumption. The concept of effort in motor control is not uniquely defined (see for example [START_REF] Berkeley | Physical principles for economies of skilled movements[END_REF] for some possi-ble acceptations). In this paper "effort" is evaluated in terms of global change in motor control variables along the sequence. GEPPETO implements this assumption with a cost function favoring small variations of the motor control variables across the planned sequence. The present section formulates a Bayesian model that aims at performing the same planning task as GEPPETO. In this Bayesian version, the cost function implementing the "minimum effort" assumption is cast as an additional constraint on the transitions between motor control variables.

The model follows the Bayesian Programming approach illustrated in the previous section. As it will be seen, the previous single-phoneme model appears to be nested as a substructure of the three-phoneme model.

Description

Variables Planning a sequence of phonemes involves the same M , S and Φ variables considered in the previous single-phoneme model. They correspond to the motor control variables, spectral properties of the acoustic signal and phonemes, respectively. However, as we are considering a sequence instead of a single phoneme, each variable is repeated as many times as there are elements in the sequence. We therefore distinguish each different instance of the variables by an index specifying its position in the sequence. Thus, variables become M i , S i and Φ i , with i ∈ {1 : 3}. For simplicity, we will denote by Y 1:3 = {Y 1 , Y 2 , Y 3 } the conjunction of different instances of a given variable Y at different positions in the sequence. An additional variable is also introduced in order to take into account the "minimum effort" assumption made in GEPPETO. This variable is denoted by C m , standing for "motor constraint". C m is a binary variable that acts as a switch, being either in the position L for "Lazy" (corresponding to the minimum effort requirement) or H for "Hyperactive" (corresponding to its opposite, a "maximum effort" requirement).

Decomposition The joint probability distribution is P (M 1:3 S 1:3 Φ 1:3 C m ). Defining X i = {M i , S i , Φ i } as the set of all the variables at position i in the sequence, the joint probability distribution can be written

P (M 1:3 S 1:3 Φ 1:3 C m ) = P (X 1 X 2 X 3 C m ).
(9) Applying the chain rule, the right term of Equation ( 9) can be decomposed as

P (X 1 X 2 X 3 C m ) = P (X 1 ) P (X 2 | X 1 ) P (X 3 |X 2 X 1 ) P (C m | X 3 X 2 X 1 ). ( 10
)
This expression can now be simplified thanks to the hypotheses made in GEPPETO. First, besides the cost function implementing the minimum effort assumption in GEPPETO, there is nothing creating any dependencies relating variables at different positions in the sequence. In real speech production, this type of constraint does exist, and corresponds to what is called "phonotactic" rules in linguistics. These rules are language dependent. It is not the purpose of the present study to address this type of high level linguistic constraints. We see though that the Bayesian Programming framework would be appropriate to account for this kind of constraint. According to this independence of variables at different positions in the sequence, the second and third factors in the decomposition of Equation (10) can be simplified such that:

P (X 1:3 C m ) = P (X 1 ) P (X 2 ) P (X 3 ) P (C m | X 3 X 2 X 1 ). ( 11
)
The last factor in the decomposition corresponds to the dependence of the variable C m on the other variables. According to H 6a , the cost function only takes into account control variables M 1:3 at all positions in the sequence by computing the perimeter of the triangle defined by these control variables. No other variable directly influences variable C m . The last term in the decomposition of Equation ( 11) can therefore be further simplified as

P (C m | X 3 X 2 X 1 ) = P (C m | M 3 M 2 M 1 ), (12)
Taking into account these simplifications, the joint probability distribution becomes P (X 1:3 C m ) = P (X 1 ) P (X 2 ) P (X 3 )

P (C m | M 3 M 2 M 1 ). ( 13
)
From this last expression it can be seen that the decomposition of the joint probability distribution P (X 1:3 C m ) contains three copies of the probability distribution P (X i ). They correspond to the joint probability distribution of variables involved in the production of a single phoneme, derived in Section 3.1:

P (X i ) = P (M i ) P (S i | M i ) P (Φ i | S i ). ( 14
)
Combining Equation ( 13) with Equation ( 14) gives the complete decomposition of the joint probability distribution:

P (M 1:3 S 1:3 Φ 1:3 C m ) = P (M 1 ) P (S 1 | M 1 ) P (Φ 1 | S 1 ) P (M 2 ) P (S 2 | M 2 ) P (Φ 2 | S 2 ) P (M 3 ) P (S 3 | M 3 ) P (Φ 3 | S 3 ) P (C m | M 3 M 2 M 1 ). ( 15
)
Figure 8 represents the Bayesian network corresponding to the decomposition of the generative model given by Equation ( 15). Although this decomposition does not show explicit dependencies between control variables, the posterior distributions -of the most likely motor sequence, given a sequence of phonemes -will be conditionally dependent. In other words, the most likely motor change towards the next phoneme depends on all the other phonemes of the sequence. This decomposition can be interpreted as being composed of the likelihoods of producing the desired phonemes at each target point in the sequence ( i.e. all P (S i | M i ) P (Φ i | S i ) terms in Equation ( 15)) and of a prior belief about the sequence of motor goals ( other terms in Equation ( 15)). We will see below that the cost function in the control space, described in H 6 , plays the role of this prior belief, while the perceptual constraints can be regarded as the corresponding likelihoods.

Parametric forms Having derived the decomposition of the joint probability distribution, Equation ( 15), it is necessary to determine the form taken by each of the factors in this expression. This was already done in

M 1 S 1 Ф 1 M 2 S 2 Ф 2 M 3 S 3 Ф 3 C m
Fig. 8 Bayesian network corresponding to the decomposition of the joint probability distribution given by Equation (15).

Section 3 for the terms appearing in the first three lines in Equation ( 15). The last term, P (C m |M 3 M 2 M 1 ), represents the dependence of variable C m on the control variables. The aim of the cost function in GEP-PETO is to penalize patterns of control variables that are far from each other by attributing them a cost that increases with the perimeter of the triangle that they define in the control space (H 6a ). The same motor constraint is implemented in

P (C m | M 3 M 2 M 1 ) through: P ([C m = L] | M 3 M 2 M 1 ) = e -κ M (|M 2 -M 1 |+|M 2 -M 3 |+|M 3 -M 1 |) . (16)
The additional parameter κ M is introduced in order to modulate the strength of the constraint on motor control variables M . The motor constraint given by Equation ( 16) is interpreted in the following way. The further the control variables are from each other, the smaller the probability for the variable C m to be in state L ="Lazy". Therefore if the state of being lazy is desired, its realization would become more probable for motor control variables being close from each other.

For completeness, as C m takes only two values, the corresponding expression for the probability of having C m = H is given by:

P ([C m = H] | M 3 M 2 M 1 ) = 1 -P ([C m = L] | M 3 M 2 M 1 ) = 1 -e -κ M (|M 2 -M 1 |+|M 2 -M 3 |+|M 3 -M 1 |) .
(17)

Planning in the context of the Bayesian three-phoneme model

Considering the planning problem addressed in GEP-PETO, the task assigned to the Bayesian three-phoneme model is to infer a sequence of motor control variables M 1:3 under the condition that the desired phonemic categories Φ 1:3 are reached and assuming the "Lazy" state for variable C m . This inference is formulated in Bayesian terms by P (M 1:3 |Φ 1:3 [C m = L]). This is again solved in a standard way through the knowledge provided by the joint probability distribution of Equation (15). The corresponding expression is given by:

P (M 1:3 | Φ 1:3 [C m = L]) ∝ S 1:3 P (X 1 )P (X 2 )P (X 3 )P ([C m = L]|M 3 M 2 M 1 ) ∝ P ([C m = L] | M 3 M 2 M 1 ) S 1:3 i=1:3 P (X i ) ∝ P ([C m = L] | M 3 M 2 M 1 ) i=1:3 S i P (M i S i Φ i ) ∝ P ([C m = L] | M 3 M 2 M 1 ) i=1:3 P (Φ i | S * (M i )) (18)
where the proportionality symbols account for normalization constants. The last line derives from the same observation as in Section 3 for the sum over S i of the joint probability distribution P (M i S i Φ i ). This completely specifies the solution to the inference problem.

Results

The Bayesian three-phoneme model is implemented through Monte Carlo sampling as described in Section 3.1.3 The performance of the previous single-phoneme model was evaluated in relation to its capacity to produce spectral properties that are located in the desired target regions in acoustic space . For the present threephoneme model, performance is evaluated with respect to the minimization of the distance between inferred control variables along with the same perceptual constraint concerning the spectral properties. With respect to this perceptual constraint, Figure 9 illustrates the spectral properties of acoustic signals corresponding to 100 inferred motor control variables for the production of the sequence /aki/. The optimal solution obtained under similar conditions with GEPPETO is also displayed for comparison.

The first observation is the variability of the results obtained from the Bayesian three-phoneme model. This was expected given the probabilistic framework of the model. Secondly, it can be observed that the obtained spectral patterns effectively distribute inside the correct target regions. This illustrates that the three-phoneme model correctly satisfies the perceptual constraint. It will be shown below that the achievement of this constraint can be controlled by the values of parameters κ S and κ M . Note that these spectral and motor parameters are controlling the certainty or confidence of probabilistic mappings and are thus related to precision (or inverse variance) of the control in the corresponding space. Thirdly, it can be noted that point clouds characterizing the distributions of the resulting spectral properties are shifted from the center of the target regions toward their boundaries, with a clear tendency for the /a/ productions to be shifted to smaller F 1 values, and for the /k/ productions to be shifted toward higher F 2 values. This shows the influence of the motor constraint on the planned sequence at the acoustic level. This is also observed in the sounds obtained with GEPPETO. Finally, Figure 10 illustrates the role of parameters κ M and κ S in the fulfillment of the perceptual constraint. It can be seen that the stronger the weight of the motor constraint (κ M ), relative to the perceptual constraint (κ S ), the stronger the shift of the points from the central regions. At the extreme, targets are no longer reached if the value of κ M becomes too large compared to κ S as can bee seen in the two bottom panels of Figure 10.

We have seen the effect of the motor constraint on the planned sequence at the acoustic level: acoustic signals deviate from the center of the target regions and tend to be closer from each other. However, it should be noted that the minimization of the motor cost occurs in the motor space and not in the acoustic space. Hence, the closer proximity of spectral realizations of the phonemes in the sequence is a consequence in the acoustic space of the constraint in the motor space. This explain in particular that, in the upper panel of Figure 10, the spectral characteristics of the selected realizations of phoneme /i/ appear to deviate away from the two other phonemes, instead of being closer as one would expect from the form of the motor constraint. A tentative explanation for this phenomenon is the strong non-linearity of the mapping relating motor control variables to the spectral properties of the acoustic signal, observed in particular for vowel /i/.

In order to evaluate whether the motor constraint actually performs the minimization of the distance between motor control variables involved in the sequence, it is necessary to evaluate the actual perimeter of the triangle that they define in the motor space. Figure 11 shows the average value taken by this perimeter for 100 inferences of the Bayesian three-phoneme model for the sequence /aki/, as a function of the parameter κ M and for different values of κ S . The value obtained with GEPPETO is also presented for comparison. We first note that curves corresponding to different values of κ S all merge for κ M = 0. This corresponds to the situation where there is no constraint in the control variables, and therefore planning of the sequence is performed independently of the other phonemes in the sequence. The average value of the distance between control variables does not depend on κ S in that case. Next, we observe that the average perimeter is clearly reduced when the strength of the motor constraint is raised with κ M , and the capacity to minimize the motor cost is stronger for higher values of κ S (i.e. for small perceptual constraints, see Figure 5). This illustrates the trade-off between the two constraints governed by κ M and κ S in the Bayesian three-phoneme model. Now, how does the Bayesian three-phoneme model perform compared to GEPPETO? It can be noted in Figure 11 that for each value of κ S (i.e. each level of perceptual constraint) there is a value of κ M (i.e. a strength of the motor constraint) for which the average distance between control variables obtained with the Bayesian three-phoneme model coincides with the result obtained with GEPPETO. For instance, for κ M = 1 the Bayesian three-phoneme model coincides with GEPPETO when κ S = 0.2. Figure 9 confirms that for these specific parameter values, the perceptual constraint is correctly satisfied and the spectral characteristics obtained with the Bayesian three-phoneme model are close to the spectral characteristics obtained with GEPPETO. This suggests the equivalence of the two models for these specific values of the parameters. However, if we compare the optimal control variables obtained with the Bayesian three-phoneme model, i.e. those that minimize the perimeter in the motor control space, with the optimal commands obtained with GEPPETO, we realize that the optimal perimeter obtained by the Bayesian three-phoneme model is actually smaller than the one obtained by GEPPETO. This suggests that GEPPETO has not found the true optimal values. We will return to this issue in the next section. Keeping κ S to the same value as in Figure 9 and augmenting κ M by a factor 5. Middle: Keeping κ M to the same value as in Figure 9 and multiplying κ S by a factor 5 (remember that augmenting κ S corresponds to relaxing the constraint, see Figure 5). Bottom: Augmenting the motor constraint and relaxing the perceptual constraint at the same time. Phonemic targets are attained as long as there is a correct balance between the strength of the motor constraint and the strength of the perceptual constraint.

Discussion

Equivalence of models

We have described a feedforward optimal control model of speech planning formulated within a Bayesian model- ing framework. The results of simulations indicate that, as for its optimal control version, the Bayesian threephoneme model correctly infers motor control variables that perform the desired motor task satisfying the specified perceptual and motor constraints. Furthermore, for specific values of the parameters characterizing the strengths of the constraints in the Bayesian three-phoneme model, simulations suggest the equivalence of results obtained by both models. This equivalence is evaluated on the basis of the comparison of average values obtained with the Bayesian three-phoneme model with the optimal solution obtained with GEPPETO. Nevertheless, it can be shown that the optimal control model can be obtained as a particular case of the Bayesian three-phoneme model if one looks for the configuration of control variables that maximize the posterior probability given by P (M 1:3 | Φ 1:3 ). The derivation of this result is provided as a supplementary material and rests on the property that the negative logarithm of P (M 1:3 | Φ 1:3 ) turns out to be equivalent to the cost function of GEPPETO. Therefore, maximizing the probability P (M 1:3 | Φ 1:3 ) is identical to minimizing the equivalent cost function of GEPPETO, showing that the Bayesian three-phoneme model can be simplified to GEPPETO in this specific implementation scheme. Note that there are mathematical theorems showing that a Bayesian scheme exists for any set of cost functions and optimal behaviour. These are known as complete class theorems [START_REF] Brown | A complete class theorem for statistical problems with finite sample spaces[END_REF][START_REF] Robert | The Bayesian Choice -From Decision-Theoretic Foundations to Computational Implementation[END_REF]. Know-ing this theoretical context, stating that the Bayesian reformulation of GEPPETO is able to account for its optimal control scheme is not surprising. However, the theorems state the existence of the Bayesian reformulation; our contribution goes further, by defining the structure taken by this reformulation in our case. This is discussed in more detail in the sections below. Note also that parameters κ S and κ M are absent from the GEPPETO model. The inference perspective on motor control equips models with parameters that encode confidence or precision. In other contexts, these parameters could reflect important sources of inter-subject variability; and, possibly, an explanation for neurological and psychiatric symptoms (e.g. Parkinson's disease). In addition, such parameters could have a key role in modulating the gain of policy selection or motor execution and may play a pivotal role in phenomena like sensory attenuation and action observation [START_REF] Friston | Action understanding and active inference[END_REF].

Addressing redundancy and variability in formal terms

We were interested in the problem of how a feedforward model of motor planning can solve the indeterminacy characterizing the specification of motor control variables for achieving a desired motor task, without resulting in a stereotyped behavior. The essence of the dilemma was rooted in the fact that on the one hand indeterminacy arises from redundancy, i.e. from the multiplicity of solutions to the problem, and on the other hand solving redundancy, i.e. eliminating all possible solutions but one, inevitably results in stereotypy. We suggested that variability could be recovered at this point by assuming that even if the planning problem is driven by an optimality assumption, the actual solution might not be a stereotyped one. The absence of stereotypy may be first due to inherent computational limitations of the search for optimal solutions. In GEP-PETO the optimization algorithm relies on a gradient descent scheme. Crucially, due to non-linearities relating variables in the model, the cost function may feature multiple local minima and the solutions obtained by gradient descent techniques may be highly dependent on the initial values of the optimization algorithm. Initializing the gradient descent algorithm in GEPPETO with different starting positions does indeed drive convergence to different locally "optimal" solutions. In particular this explains why the solution obtained with GEPPETO, as shown in Figure 9, appears to have a greater perimeter value than the optimal solution found with the Bayesian three-phoneme model. The result for GEPPETO was actually chosen as the best one out of 100 different initializations of the descent algorithm.

The fact that the gradient descent algorithm has failed to converge in all of these 100 initializations indicates the degree of complexity of the optimization process.

In this context, it could be argued that variability in speech production arises from the existence of these multiple local solutions into which the optimization process may differently converge depending on its initial configuration. However, the variability introduced this way cannot be formally justified as actually arising from the model itself, since it is just an indirect consequence of the failure of its implementation for finding the true optimal solution that the model actually predicts. Moreover, this ad hoc implementation can only account for variability in a qualitative way and does not have any theoretical or cognitive foundation.

In contrast, formulating the feedforward planning process within a Bayesian modeling framework has allowed us addressing the indeterminacy of the problem in addition of dealing with variability in formal terms. This is made possible by the fact that the Bayesian approach does not solve indeterminacy by suppressing all solutions but one. Instead, the Bayesian framework characterizes every possible configuration by its probability to achieve the task. Redundancy is then solved by randomly selecting motor control variables under the corresponding probability distribution. The optimal achievement of the task is still ensured in average, since the most probable motor control variables inferred under this process correspond to the more relevant ones for the motor task. Variability becomes therefore an inherent consequence of the formalism. Furthermore, the variability generated with this approach has a specific structure that could be compared with experimental data. For instance the model predicts that the relative frequency of selected motor control variables is given by the probability P (M 1:3 | Φ 1:3 ).

Therefore, the advantage of the Bayesian modeling approach is to suggest that a probabilistic description of the planning process is able to deal with the selection of solutions to an ill-posed problem without destroying variability [START_REF] Brown | A complete class theorem for statistical problems with finite sample spaces[END_REF]. This allows to treat variability in formal terms and not as the result of an ad hoc implementation of the model.

The pertinence of an approach that designs models integrating multiple local solutions in formal terms is illustrated by the work of [START_REF] Ganesh | Motor memory and local minimization of error and effort, not global optimization, determine motor behavior[END_REF]. Their work indicates that motor memory plays a crucial role influencing the outcome of the planning process, in addition to the optimization of cost related to error and effort. Thus, motor memory would be responsible for setting variable initial states of the motor system, which would influence the convergence of the search for opti-mal solutions toward local optima. Even if the Bayesian three-phoneme model that we have presented does not account for the role of motor memory in the planning process, the Bayesian modeling approach offers a framework in which motor memory could be modeled via a set of local approximations to the complete probability distribution, as it would be performed by local Laplace approximation or by standard variational inference methods. This raises the question of how agents would encode the knowledge described by the probability distributions involved in the presented scheme. While a complete representation of a complex knowledge would involve an important amount of resources, it would be natural to select a simpler approximation to this knowledge as it would be advantageous for the agent and often sufficient for practical purpose. Indeed, there is a fairly established literature on active inference using variational Bayes in the context of Bayesian filtering (also known as predictive coding). In brief, by equipping predictive coding with reflexes, active inference simulates action trajectories, action observation and indeed communication (e.g., the bird song examples in [START_REF] Friston | Active inference, communication and hermeneutics[END_REF]). As in the current Bayesian formulation, active inference dispenses with cost functions and replaces them with prior beliefs about the way the motor plant should behave.

Conclusion

We propose to conclude by widening the discussion of our contribution. We first note that, of course, marrying optimal control or optimal planning theories and probabilistic modeling has already a long history. Previous approaches abound; we provide a few, mostly classical entry points in the vast literature of decisiontheoretic planning in robotics and AI [START_REF] Kaelbling | Planning and acting in partially observable stochastic domains[END_REF][START_REF] Boutilier | Decision theoretic planning: Structural assumptions and computational leverage[END_REF][START_REF] Attias | Planning by probabilistic inference[END_REF]Murphy, 2002;[START_REF] Toussaint | Probabilistic inference as a model of planned behavior[END_REF][START_REF] Kappen | Optimal control as a graphical model inference problem[END_REF], and in the Bayesian Decision Theory in cognitive modeling [START_REF] Wolpert | Probabilistic models in human sensorimotor control[END_REF][START_REF] Daunizeau | Observing the observer (I): Meta-bayesian models of learning and decision-making[END_REF][START_REF] Ma | Signal detection theory, uncertainty, and poisson-like population codes[END_REF][START_REF] Ma | Organizing probabilistic models of perception[END_REF] However, our approach differs in that, instead of marrying probabilities and cost functions, we proposed to do away completely with the notion of cost functions , something that has already been proposed within the active inference scheme [START_REF] Friston | What is optimal about motor control?[END_REF][START_REF] Friston | Reinforcement learning or active inference?[END_REF][START_REF] Friston | Active inference and agency: optimal control without cost functions[END_REF]. We have reformulated the cost function of an existing, optimality based model, as probability distributions, motivated by a desire to obtain trial-to-trial variability in a principled manner. A consequence that we have already exposed is that the optimality based model could be seen as a special case of the proposed Bayesian three-phoneme model. However, this nesting of models does not imply that the Bayesian modeling framework would be inherently more powerful than optimal based modeling. Indeed, one could strive to expand optimal based models to recover trial-to-trial variability, with other mathematical methods that we do not imagine at the moment. Therefore, this warrants caution concerning the interpretation of our contribution. Since a single model can be expressed equivalently in two different mathematical formalisms, none of these formalisms can be claimed, at face value, to be e.g. more biologically plausible than the other. This is a singular epistemological stand, in the current debate about the theoretical contribution of Bayesian modeling to cognitive sciences [START_REF] Jones | Bayesian fundamentalism or enlightenment? on the explanatory status and theoretical contributions of bayesian models of cognition[END_REF][START_REF] Bowers | Bayesian just-so stories in psychology and neuroscience[END_REF][START_REF] Hahn | The Bayesian boom: good thing or bad?[END_REF]. A single Bayesian model does not bring much evidence that the brain would encode and manipulate probabilities (the so-called Bayesian Brain Theory), if only because of mathematical equivalent model found in other formalisms. Indeed, the Bayesian formalism itself appears under-constrained. In other words, writing a Bayesian model of a given cognitive function is always feasible. As a side note, this does not preclude a reifiability based argument in favor of the Bayesian Brain Theory; if there are many Bayesian models of many cognitive functions, then one can find probabilities more likely to be "used" by natural cognitive systems. Our current contribution is a step in this research program. This is also complementary to, and orthogonal to, studies of the same question at the microscopic biological level (e.g., neural or population of neuron based accounts of cognitive processes; [START_REF] Pouget | Probabilistic brains: knowns and unknowns[END_REF]).

Advantages, then, are to be found on other grounds. We propose to highlight the interest of Bayesian modeling as a mathematical modeling tool. As we have seen, to express knowledge, a single mathematical construct, that is, probability distributions, is required. Such a unified formalism is interesting in several aspects. For instance, it makes comparison and composition of pieces of knowledge easier. In our case, this was illustrated by the composition of speech constraints of varied nature, one concerning motor economy, the other concerning perceptual discriminability. Furthermore, additional constraints are easy to combine in order to enrich the model (e.g., we added to the model, but did not describe here for brevity, a parallel branch to acoustic targets that controls the force output, in order to modulate speech rate while conserving intelligibility).

Interpretability of the model also benefits from this unified formalism. In the present Bayesian model, the perceptual step-shaped constraint of GEPPETO was derived from the inversion of an internal representation of phonemes in terms of simple distributions of the spectral properties of the produced sounds. The perceptual constraint is therefore interpreted in the Bayesian model as an internal perceptual categorization of the produced acoustic sound. This is interesting as it shows explicitly how perception is assumed to be involved in the control process. It further illustrates an additional advantage of the Bayesian framework as being well suited for treating perception and action in a unified framework [START_REF] Toussaint | Probabilistic inference as a model of planned behavior[END_REF][START_REF] Friston | The free-energy principle: a unified brain theory?[END_REF].

Moreover, to manipulate knowledge, the few rules of probability calculus are sufficient. Inference directly and automatically derives from the choice of the model; in more technically precise words, defining the joint probability distribution and the probabilistic questions asked to this joint probability distribution completely constrains the resulting inference processes. In that sense, our method differs from most other approaches, by placing the focus of modeling on knowledge expression, instead of the inference process. The cognitive model we propose is therefore resolutely representational; it lies at the algorithmic level of Marr's classical hierarchy [START_REF] Marr | Vision. A Computational Investigation into the Human Representation and Processing of Visual Information[END_REF]. We expect this original perspective to expand in cognitive science modeling, in the wake of the current explosion of probabilistic programming languages [START_REF] Goodman | Church: a language for generative models[END_REF][START_REF] Gordon | Probabilistic programming[END_REF].

Fig. 1 Fig. 2

 12 Fig. 1 Projections of the 3-dimensional dispersion ellipsoids corresponding to each target region characterizing phonemes. Top: (F 2 , F 1 ) plane; Bottom: (F 2 , F 3 ) plane. The dotted line on the bottom image indicates the F 3 value specified in Figures 4 and 5. (/oe/ is noted /A/ in all figures for sake of simplicity in graphic representation)

Fig. 3

 3 Fig.3Bayesian network representing the decomposition of the joint probability distribution given by Equation (4).

Fig. 4

 4 Fig. 4 Values of the probability distribution P (S | [Φ = k]) (left) and the corresponding likelihood function P ([Φ = k] | S) (right), projected on the (F 1 , F 2 ) plane defined by F 3 = 2, 450 Hz.

Fig. 5

 5 Fig. 5 Effect of the κ S parameter on the likelihood functions P (Φ | S). We superpose the likelihood functions for all phonemes, and project them on the plane (F 1 , F 2 ) defined by F 3 = 2, 450 Hz. Top: κ S = 0.3. Bottom: κ S = 1. Right panels are top-views of the left panels. Smaller values for κ S narrow the confidence regions in the categorization task.

λ 3

 3 in phoneme /O/. Some other appear to have a wide range of variation, for instance λ 5 for all phonemes. This indicates the importance of the role of each muscle in performing each phoneme. In particular, we notice that control variables λ 1 and λ 3 negatively correlate for phonemes /i, e, E, a/. Smaller values of λ 1 are related to higher values for λ 3 and vice versa. The values taken by these control variables specify the activation level of the Posterior Genioglossus and Hyoglossus muscles.

Fig. 7

 7 Fig. 7 Histograms of the acoustic signals resulting from the sampling of control variables sampled from the inferred probability distribution P (M | Φ) for κ S = 0.2. The vertical dotted lines indicate borders of the acoustic regions characterizing each phoneme in the original formulation of GEPPETO.

F2Fig. 9

 9 Fig.9Projection of the acoustic signal on the (F 2 , F 1 ) plane, obtained by 100 motor control variables sampled from the inference probability distribution for the production of sequence /aki/. The acoustic signal obtained with GEPPETO and by the sample of the Bayesian three-phoneme model with minimum perimeter are also indicated. Values of the perimeters obtained by each model are indicated.

F2Fig. 10

 10 Fig. 10 Effect of parameters κ S and κ M on the spectral properties of the acoustic signals obtained by the Bayesian three-phoneme model. Refer to Figure 9 for comparison. Top:Keeping κ S to the same value as in Figure9and augmenting κ M by a factor 5. Middle: Keeping κ M to the same value as in Figure9and multiplying κ S by a factor 5 (remember that augmenting κ S corresponds to relaxing the constraint, see Figure5). Bottom: Augmenting the motor constraint and relaxing the perceptual constraint at the same time. Phonemic targets are attained as long as there is a correct balance between the strength of the motor constraint and the strength of the perceptual constraint.

Fig. 11

 11 Fig. 11 Average distances obtained with the Bayesian three-phoneme model as a function of the parameter κ M . Results for different values of κ S are plotted (listed in the insert on the right). Error bars indicate variability obtained over 100 random samplings. The black horizontal dashed line represents the value obtained with GEPPETO.

Acknowledgements Authors wish to thank Pierre Bessière and Jean-Luc Schwartz for guidance and inspiring conversations.

The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 339152, "Speech Unit(e)s", PI: Jean-Luc-Schwartz).


QUASI REAL-TIME SUB-SPACE 3D DEFORMABLE FUSION

Bharat Singh, Stavros Alchatzidis and Nikos Paragios

Ecole Centrale Paris

ABSTRACT

In this paper we propose a novel and robust approach for deformable fusion using a metric defined in an appropriate subspace which is adaptive to the image-content/image-modality. We adopt a graph-based formulation that assumes that intensities of corresponding pixels in the two image domains are related through an unknown piece-wise constant linear function. This relation is propagated to an appropriate sub-space (wavelets coefficients) where a criterion that couples the estimation of the deformation field with optimal transport function on the subspace and the smoothness of the deformation is considered. Message passing methods efficiently implemented using jump flooding are considered to get the optimal parameters of the deformation field and the transport function. Promising results are obtained in near real-time with a unique parameter setting.

Index Termsdeformable image registration, wavelets, similarity metric

INTRODUCTION

Deformable multi-modal fusion is a well explored topic in biomedical image analysis. It shares the challenges of mono-modal registration with the additional constraint of defining an appropriate similarity metric capturing the codependencies between the images to be registered.

Removing intensity variation through implicit image normalization was the first attempt to cope with multi-modal fusion using criteria like normalized cross correlation or correlation ratio that work well in general in the multi-modal setting. Exploring feature spaces beyond intensities that capture rich geometric information was a further attempt to cope with deformable fusion [START_REF] Ou | DRAMMS: Deformable registration via attribute matching and mutual-saliency weighting[END_REF]. The use of statistical methods and in particular of mutual information [START_REF] Viola | Alignment by maximization of mutual information[END_REF] was a major breakthrough in the field because, the similarity criterion does not measure the quality of correspondences between individual pixels but a statistical entropy-driven criterion derived from the distributions of the image intensities before and after alignment. Numerous metrics have been inspired from mutual information, for example its normalized version, the Kullback-Leibler divergence or the Jensen-Renyi divergence. It should be noted that due to the non-linearity of the similarity criterion, the optimization component of the registration process often converges to a local minimum and the obtained solution depends heavily on the initial conditions/alignment between the two images.

The aim of this paper is to provide a novel, efficient, robust metric for general purpose multi-modal registration. To this end, we consider a robust/efficient sub-space local representation of the images on a wavelet base. Then assuming a local partition of the image domain, where there exists an unknown piece-wise constant linear mapping from one modality to another, we propagate this mapping to the subspace. An appropriate metric that depends on the unknown mapping between the two images is defined in the subspace. Such a propagation makes feasible the quantization of the mapping parameter through the subspace metric, which is propagated to a hierarchical/coarse-to-fine graphical model [START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF]. In this model both deformations (grid-like interpolations [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF]) and local mapping are unknown and optimized together. The resulting MRF is sub-modular and message passing methods are considered to recover the joint solution of subspacemapping/deformation producing quasi realtime results which are comparable to the state of the art.

The reminder of this paper is organized as follows, section 2 presents the metric and the registration formulation, while section 3 describes the the optimization framework. Implementation details and experimental validation are part of the section 4 while section 5 concludes the paper.

DEFORMABLE FUSION

Let us consider two images s(x) and t(x) corresponding to different modalities. Without loss of generality let us consider that the aim is to map s to t using a non linear transformation d(x). Once a similarity metric is given, denoted with ρ(s(x), t(y)), the optimal solution of fusion refers to the transformation d(x) that minimizes an image-driven term and a regularization constraint on the deformation space,

min d Ω ρ(s(d(x)), t(x))dx + α Ω φ(∇d(x))dx
where Ω is the image domain and φ imposes a smoothness constraint on the deformation field. Let us consider a partition of the image domain Ω into a set of sub-domains

{Ω 1 , • • • Ω n }.
The union of these domains correspond to the image registration domain and we further assume that adjacent sub-domains could overlap. Then, given the optimal deformation, we consider a linear transport function relating the two signals, x ∈ Ω i : s(d(x)) = a i t(x) + b i , where a i and b i are sub-domain constants corresponding to the linear mapping between the two modalities. It should be noted that linearity is for the purpose of facilitating the introduction of the model and by no means a hard constraint. Given such a model, we can cast registration as a coupled inference process that aims at optimally determining the image partition sub-domains, the transport coefficients and the optimal deformation, or

min n,ai,bi,Ωi,d n i=1 Ωi ρ(s(d(x)), a i t(x) + b i )dx +α Ω φ(∇d(x))dx

Deformation Model/Image Partition

We consider a grid like representation endowed with an interpolation function. Let us consider the well established model of hierarchical free-form deformations which consists of superimposing a grid-like structure

G k,l , (k, l) ∈ [1, K] × [1, L]
to the image. The deformation of a given image pixel is computed by a cubic B-Spline interpolation of the displacements of the "adjacent" grid nodes [START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF].

The advantage of such a deformation is that it significantly reduces the degrees of freedom in registration, guarantees smoothness of the deformation field and provides a natural partition of the image domain. Let us consider that each control point is endowed with a rectangular support domain R k,l that is centered to the control point. In such a case the registration functional can be expressed as

min a k,l ,b k,l ,d k,l K k=1 L l=1 Ri ρ(s(d(x)), a k,l t(x) + b k,l )dx +α Ω φ(∇d(x))dx
The number of sub-domains, and their corresponding integrals are not considered as unknowns.

Similarity Metric

Without loss of generality, let us consider a patch around a control point in a 1-D signal with a support domain R k . Let us also assume that the control point is in the source image and its intensity distribution is given by s(d(x)). The intensity distribution of the corresponding patch in the target image is given by a k t(x) + b k . Our aim is to define a similarity criterion which is minimized when the source image aligns with the target image, irrespective of a k and b k . In a discrete domain, let us assume the length of the source and the target to be N . For a signal f of size N , the discrete Haar wavelet transform is defined as,

dwt[i] = f [2i]+f [2i+1] √ 2 , when 0 ≤ i < N 2 and dwt[i] = f [2i-N ]-f [2i+1-N ] √ 2 , when N 2 ≤ i < N
, where dwt represents the first level wavelet coefficients. Such an operation is applied recursively on the first half of the signal. It is clear that all the coefficients (other than the first level wavelet coefficients, LL1 in case of 2 levels) comprise of a local difference of intensities in the input signal at various scales. In the process, the constant a k is scaled up and b k gets canceled. Applying discrete wavelet transform over the support domain of the control point, would lead to wavelet coefficients of the form s (d(x)) and a k t (x), where s and t are the wavelet coefficients of the source and the target respectively. A linear relationship between the intensities of the source and the target images would also be propagated to the wavelet space. The relation between wavelet coefficients would then be given by, s (d(x)) = a k t (x). After taking magnitude of the coefficients and applying logarithm, we have log(|s

(d(x))|) = log(|t (x)|) + log(a k ).
Using this property, we define the similarity metric ρ(s(d(x)), t(x))) as,

1 N N i=1 (log(|s (d(x))|) -log(|t (x)|) -log(a k ))
In the metric which we have defined, the factor due to linear intensity variation appears as a constant instead of a multiplicative factor. This constant is small when compared to scale of the signals and is spatially stationary. Therefore, it does not alter the nature of the distribution of the similarity metric. Towards accounting for more complex signal variations, we assume a non-linear exponential relationship between the wavelet coefficients, s (d(x)) γ k = a k t (x), where γ k is the exponential variation in the wavelet coefficients. To compensate for this non-linear change, we define the similarity cost ρ (s(d(x)), t(x)), c k ) as,

1 N N i=1 (c k • γ k • log(|s (d(x))|) -log(|t (x)|) -log(a k ))
where ρ is the joint distribution of the similarity metric with the non-linear exponential constant c k .

GRAPH BASED MULTI-MODAL SUBSPACE FUSION

Inspired by the work of [START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF], we reformulate registration as a discrete labeling problem. Let us consider a discrete set of labels L = {l 1 , l 2 , ..., l i } which correspond to a quantized version of the deformation space θ = {d 1 , ..., d i } and a second label set

⊗ = {ω 1 , • • • , ω j } corresponding to a quantized version of the non-linear transport coefficient c k , Γ = {γ 1 , • • • , γ j }.
Based on the label assigned to each In this context, registration can be formulated as a labelling problem which can be efficiently solved using markov random fields [START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF] as

E M RF (l, ω) = p∈G V p (l p , ω p ) + p∈G q∈N (G)
V pq (l p , l q )

The unary potential corresponding to the deformation and the non-linear correction constant c k is given by

V p (l p , ω p ) = Ωp ρ (s(d lp + x), t(x), γ ωp )dx
In order to impose smoothness on the deformation field, we use the pair-wise terms of the graphical model as mentioned in [START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF]. For imposing local-consistency on the non-linearity correction, we can impose on the pair-wise term a Potts like model between the distance of the retained coefficients. However, due to the robustness of the metric and its hierarchical nature such a constraint is not necessary. For recovering the optimal parameters of the MRF, message passing methods efficiently implemented using jump flooding [START_REF] Alchatzidis | Efficient parallel message computation for map inference[END_REF] are considered. Further, we use a limited label set and perform several optimization cycles. The control grid is reset after every optimization cycle and the spacing between different labels is reduced by a factor to capture finer displacements.

EXPERIMENTS AND RESULTS

The validation of deformable multi-modal registration is a challenging task since the the deformation field for the images is generally not available. We demonstrate the performance and the robustness of the proposed method by providing statistical and visual results on mono-modal and multimodal datasets.

Inter Subject Deformable Brain Registration

We performed 306 pairwise registrations on 18 MR-T1 brain images available at Internet Brain Segmentation Repository (http://www.cma.mgh.harvard.edu/ibsr/). In order to evaluate performance, we deform the segmentation available for 3 different brain parts, and compute the Jaccard Overlap between the deformed source segmentation and the target segmentation. From Fig. 1 we can infer that, our metric performs better in terms of median scores among all metrics for "Cerebrospinal Fluid" and "Gray Matter", while being at par with other metrics on "White Mater". In order to compare our approach with state of the art methods, we evaluate 10 registration algorithms [START_REF] Jenkinson | A global optimisation method for robust affine registration of brain images[END_REF][START_REF] Andersson | FNIRT-FMRIBs non-linear image registration tool[END_REF][START_REF] Roger P Woods | Rapid automated algorithm for aligning and reslicing PET images[END_REF][START_REF] Babak A Ardekani | Quantitative comparison of algorithms for inter-subject registration of 3D volumetric brain mri scans[END_REF][START_REF] Brian B Avants | Symmetric diffeomorphic image registration with cross-correlation: Evaluating automated labeling of elderly and neurodegenerative brain[END_REF][START_REF] Rueckert | Nonrigid registration using free-form deformations: application to breast MR images[END_REF][START_REF] Vercauteren | Diffeomorphic demons: Efficient non-parametric image registration[END_REF][START_REF] Ou | DRAMMS: Deformable registration via attribute matching and mutual-saliency weighting[END_REF][START_REF] Glocker | Dense image registration through MRFs and efficient linear programming[END_REF], on LONI-LPBA40 dataset [START_REF] Shattuck | Construction of a 3D probabilistic atlas of human cortical structures[END_REF]. 1560 pairwise registrations on 40 skull stripped brain images (181x217x181) were performed and average Jaccard Overlap over 56 segmentation labels was computed. Fig. 2. It can be seen that the proposed approach is comparable to state of the art methods in terms of accuracy while being orders of magnitude faster.

Multi-Modal Brain Registration

In order to test our metric on multi-modal images, we performed registration on multi-modal brain images from different modalities like MR-T1, MR-T2, CT, PET. The images are part of the online dataset provided by Vanderbilt University (http://www.insight-journal.org/rire/). Registration results of our algorithm are shown in Fig. 2.

Implementation Details

We tested our approach on a Nvidia Geforce GTX 480 GPU. For estimation of γ, c was varied from 0.1 to 4.9 with increments of 0.15. We add 1 to the absolute value of wavelet coefficients to ensure that the logarithm is positive. In case of mono-modal images it is not necessary to estimate γ. Our metric is free from sequential dependencies and is easily parallelizable. Parallelism is restricted in information theoretic measures due to atomic operations in histogram binning. Further, the computational complexity of our metric is independent of grid resolution unlike information theoretic measures, where run time increases at a finer grid resolution (due to fixed size joint histogram computation).

Experimental validation results presented-other than the proposed method-are courtsey, Ou, Y, 2012. Development and validations of a deformable registration algorithm for medical images

DISCUSSION

In this paper we have proposed a novel intensity based similarity metric for mono/multi-modal image registration that once combined with a graphical model and efficient discrete optimization is comparable to the state of the art in terms of accuracy while being an order of magnitude faster. The experimental validation shows that the trade-off between the performance/complexity of the proposed method is exceptional. The underlying metric computation complexity is O(n) and benefits greatly from implementation on modern parallel architectures (GPU). Experimental results on representative multi-modal cases demonstrate the interest in our approach. The exploration of alternative sub-space representations as well as more complex transport functions can be an interesting theoretical extension of the method.

Fig. 1 .

 1 Fig. 1. Jaccard Overlap for different metrics on IBSR dataset

Fig. 2 .

 2 Fig. 2. Average Jaccard Overlap over 56 ROIs for different registration algorithms is shown in the figure on the left. We show the running time of different algorithms in center figure. The red dots indicate the performance of the proposed method. Qualitative registration results on multi-modal images are shown in the figure to the right.


E Destyl
email: edes.destyl@univ-ag.fr

P S Nuiro
email: paul.nuiro@univ-ag.fr

P Poullet
email: pascal.poullet@univ-ag.fr

On the global behaviour of solutions of a coupled system of nonlinear Schrödinger equations

Keywords: Coupled nonlinear Srchödinger equations, Manakov model, Parity-time symmetry, finite element, solitons 2000 MSC

. The result is illustrated by numerical results which have been obtained with a finite element solver well adapted for that purpose.

Introduction

Hereafter, a basic model of propagation of weakly dispersive waves is considered by a coupled nonlinear Schrödinger (NLS) equations, which reads as follows:

ıu t = -u xx + κv + ıγu -(g 11 |u| 2 + g 12 |v| 2 )u, ıv t = -v xx + κu -ıγv -(g 12 |u| 2 + g 22 |v| 2 )v, (1)
with the coefficients of the nonlinear parts being real. The ı is the complex such that ı 2 = -1 and the coefficients κ and γ are positive constants that characterize gain and loss in wave components. This model is known for its pertinence for several applications of nonlinear optics and has been studied by several authors [START_REF] Driben | Stability of solitons in parity-timesymmetric couplers[END_REF][START_REF] Abdullaev | Zeno effect and switching of solitons in nonlinear couplers[END_REF][START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF].

Following the study of [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF], let us recall that the solution of the system (1), obeys a parity-time symmetry as soon as there is an equality between coefficients g 11 = g 22 . Actually, the parity symmetry is defined by the mapping: P(u, v) = (v, u) and the time reversal operator is defined by the following map:

T (u(t), v(t)) = (ū(-t), v(-t)),

where the notation ū(t) represents the complex conjugate of u(t).

Then, the nonlinearity of the system (1) is such that for any solution (u(t), v(t)) defined in a symmetric interval I = [-t 0 , t 0 ], there exists another one (u PT (t), v PT (t)) defined in I, (u PT (t), v PT (t)) = (v(-t), ū(-t)).

Main results

According to [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF], one knows the existence of a unique global solution (u(t), v(t)) ∈ C(IR, (H 1 (IR)) 2 ) of the Cauchy problem for the generalized Manakov system (1) with (u(0), v(0)) = (u 0 , v 0 ) ∈ (H 1 (IR)) 2 .

In the Hamiltonian case (for γ = 0), it is known that for (u, v) ∈ C(IR, H 1 (IR) × H 1 (IR)), the two following quantities are conserved, namely the density,

Q(t) = I R (|u| 2 + |v| 2 )dx,
and the total energy of the system (1)

E(t) = I R |u x | 2 + |v x | 2 + κ(ūv + uv) - g 11 2 |u| 4 - g 22 2 |v| 4 -g 12 |u| 2 |v| 2 dx
Our next theorem results in another technique which allows us to find an estimate of the semi-norm in H 1 of the solution of the generalized Manakov system, sharper than that can be obtained from the previous study [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF].

Theorem 1. Assuming that g 11 = g 12 = g 22 then, for the unique solution (u(t), v(t)) ∈ C(IR, (H 1 (IR)) 2 ) of the Cauchy problem for the generalized Manakov system (1) with (u(0), v(0)) = (u 0 , v 0 ) ∈ (H 1 (IR)) 2 , there exists a positive constant C such that the semi-norm of the solution has the following upper bound, for t ≤

T T 0 ( u x 2 L 2 + v x 2 L 2 )dt ≤ 4g 2 C 2 Q(0) 3 e 6γT + 2Q(0)e 2γT (κ + π T ω + 3 ) + 2πQ(0) (2)
Proof. The proof started similarly to those given in [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF], by the fact that the density and the total energy are well defined in the energy space, and that as soon as γ = 0 the integral quantities are not constant in time. First of all, by using the duality between the Sobolev spaces H -1 and H 1 , we apply the first equation of system (1) to u and the second one to v to form the density.

Then, adding the two previous dual products of left hand side of system (1) and thanks to an integration by parts, we get from the imaginary part that

d dt I R (|u| 2 + |v| 2 )dx = 2γ I R (|u| 2 -|v| 2 )dx.
That gives a control to the evolution of the density

dQ dt = 2γ( u 2 L 2 -v 2 L 2 ) ≤ 2γQ(t), (3)
and by the Grönwall inequality, the density does not blow up in finite time:

∀t ∈ [-t 0 , t 0 ], Q(t) ≤ Q(0) exp(2γ|t|). (4)
And from the real part, we obtain that

I R |u x | 2 + |v x | 2 dx + I R Re{ı(u t ū + v t v)}dx = -κ I R (vū + uv)dx + I R g(|u| 4 + |v| 4 + 2|u| 2 |v| 2 )dx. ( 5
)
Let us focus onto the second integral of the left hand side. As we assume that u can not vanish, we consider that during a time interval its imaginary part is not null (or it can be replaced by its real part). Then, using the same argument for v, we write that

Re{ı(u t ū + v t v)} = |u| 2 arctan Re u Im u t + |v| 2 arctan Re v Im v t . (6)
By integration by parts and thanks to the absolute continuity of their representants,

T 0 Re{ı(u t ū + v t v)}dt = |u| 2 arctan Re u Im u + |v| 2 arctan Re v Im v t=T t=0 - T 0 arctan Re u Im u |u| 2 t + arctan Re v Im v |v| 2 t dt.
From the absolute continuity of the representants of u and v and thanks to the oscillating behaviour of Q ( [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF]) and relation [START_REF] Driben | Stability of solitons in parity-timesymmetric couplers[END_REF], one can compute explicitly

u 2 L 2 and v 2 L 2 .
So, one obtains that they are oscillating functions also:

u 2 L 2 = 1 2γ γκC ω 2 + (γA 1 + ωA 2 ) cos 2ωt + (γA 2 -ωA 1 ) sin 2ωt , v 2 L 2 = 1 2γ γκC ω 2 + (γA 1 -ωA 2 ) cos 2ωt + (γA 2 + ωA 1 ) sin 2ωt .
Then, easy computations provide that during each time period the time derivative of u 2 L 2 and v 2 L 2 vanishes at two specific times for each component. We denote the specific instants by t

u , t (2) u = t (1) u + π/(2ω) and t (1) v , t (2)
v = t (1) v + π/(2ω). Setting α = γ/ω, and β = A 1 /A 2 , one obtains t (1) u = 1 2ω arctan αβ -1 α + β .
If α = β, one has the expression for t

(1) v , t (1) v = 1 2ω arctan αβ + 1 α -β ,
and if α = β one gets:

t (1) v = π 4ω .
An interesting reader can sort the zeros of the time derivative of u The arctan function being bounded,

I j ≤ π 2 I R t 2 t 1 |v| 2 t -|u| 2 t + t 3 t 2 |u| 2 t -|v| 2 t ≤ π 2 I R |v(t 2 )| 2 -|u(t 2 )| 2 -|v(t 1 )| 2 + |u(t 1 )| 2 +|u(t 3 )| 2 -|v(t 3 )| 2 -|u(t 2 )| 2 + |v(t 2 )| 2 dx ≤ π 2 (Q(t 1 ) + 2Q(t 2 ) + Q(t 3 )).
Obviously, the way that has been generated the upper bound for this sign configuration can be extended to any other configuration. Hence, combining and adding the upper bound for each term, one gets that,

I ≤ π 2 Q(t 0 ) + Q(t 1 ) + p-1 j=0 (Q(t 1 ) + 2Q(t 2 ) + Q(t 3 )) .
Thanks to the density estimate (4), one obtains

|I| ≤ π 2 (Q(0) + (4p + 5)Q(0)e 2γT ).
Returning to the integral of the second term of equation ( 6) left hand side, the upper bound is as follows.

T 0 I R Re{ı(u t ū + v t v)}dx dt ≤ π 2 Q(0)(1 + e 2γT ) + π 2 Q(0)(1 + (4p + 5)e 2γT ) ≤ πQ(0)(1 + T w + 3 e 2γT ).
Setting

D(T ) = T 0 ( u x 2 L 2 + v x 2 L 2
)dt, and thanks to the Gagliardo-Nirenberg inequality which is valid for function of H 1 Sobolev space [START_REF] Nirenberg | On functions of bounded mean oscillation[END_REF], we then obtain that,

D(T ) + T 0 I R Re{ı(u t ū + v t v)}dx dt ≤ κ T 0 Q(t) dt + 2gC T 0 ( u 3 L 2 u x L 2 + v 3 L 2 v x L 2 ) dt. (7)
Using the relation

a 3 b + c 3 d ≤ (a 2 + b 2 ) 3/2 (b 2 + d 2 ) 1/2
, one obtains from the inequality [START_REF] Nirenberg | On functions of bounded mean oscillation[END_REF], Jensen inequality (concavity of the square root function) and relation (4) that, D(T ) -2gCQ(0) 3/2 e 3γT D(T ) 1/2 -κQ(0)e 2γT -πQ(0)(1 + T w + 3 e 2γT ) ≤ 0, that can be considered as an univariable quadratic function of unknown D(T ) 1/2 . The discriminant is computed as

∆ = 4g 2 C 2 Q(0) 3 e 6γT + 4Q(0)e 2γT (κ + π T w + 3 ) + 4πQ(0),
which is undoubtedly positive as the density. Thus, computing the solutions, the variable D(T ) 1/2 being positive, it must stay clustered as follows:

D(T ) 1/2 ≤ gCQ(0) 3/2 e 3γT + g 2 C 2 Q(0) 3 e 6γT + Q(0)e 2γT (κ + π T ω + 3 ) + πQ(0).
Taking the square, we finally get the upper bound for t ≤ T that completes the proof.

D(T ) ≤ 4g 2 C 2 Q(0) 3 e 6γT + 2Q(0)e 2γT κ + π T ω + 3 + 2πQ(0). ( 8
)
Remark 1. The previous theorem provides an upper bound of the velocity norm in L 2 (0, T ; H 1 ) as roughly a constant times exp(6γT ). One can not compare straightly with the results that have been obtained by other studies. But, if we wish to derive from the previous study [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF], an approximation in L 2 (0, T ; H 1 ) norm, one has to combine their inequality (21) with their equality (20) and the behaviour of the mass (11). And it likely seems that this combination provides an upper bound as a constant times exp(10γT ) which is less sharp than the bound that we have obtained.

Theorem 2. Let us assume that g 11 = g 22 = g 12 .

• if γ < κ, for (u 0 , v 0 ) ∈ (H 1 (IR)) 2 a local solution at t = 0 of the system (1), there exists a constant Q max > 0 such that the solution of the (1) verifies:

sup t∈I R ||u(t)|| 2 L 2 + ||v(t)|| 2 L 2 ≤ Q max . (9)
• if γ ≥ κ, there exists a global solution of the system (1) such that,

lim t→∞ ||u(t)|| 2 L 2 + ||v(t)|| 2 L 2 = ∞.
The proof is given in [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF].

Theorem 3. Let us assume that g 11 = g 22 = g 12 .

• if γ < κ, for (u 0 , v 0 ) ∈ (H 1 (IR)) 2 a local solution at t = 0 of the system (1), there exists a constant E max > 0 such that the solution of the (1) verifies:

T 0 (||u(t)|| 2 H 1 + ||v(t)|| 2 H 1 )dt ≤ E max , (10)
• if γ ≥ κ, there exists a global solution of the system (1) such that,

lim t→∞ ||u(t)|| 2 H 1 + ||v(t)|| 2 H 1 = ∞. (11)
Proof. As g 11 = g 22 = g 12 and γ < κ and according to the conclusions of theorem 2, the density Q(t) est bounded, i.e. there exists a positive constant

Q max such that sup t∈I R ||u(t)|| 2 L 2 + ||v(t)|| 2 L 2 ≤ Q max .
Moreover the estimate [START_REF] Dias | Supercritical blowup in coupled parity-time-symmetric nonlinear Schrödinger equations[END_REF] shows that D(T ) exists. Consequently, the solution of the generalized Manakov system lies in L 2 (0, T ; H 1 ), i.e. there exists a positive constant E max such that

T 0 (||u(t)|| 2 H 1 + ||v(t)|| 2 H 1 )dt ≤ E max .
The case γ ≥ κ is a direct consequence of conclusion of theorem 2, that ends the proof.

Numerical results

The results of theorem 3 indicate that for γ < κ, the global solutions of Parity-Time symmetric system (1) remain bounded in L 2 (0, T ; H 1 ). In order to confirm this result and to investigate the behaviour of the solutions in other cases, numerical solutions have been computed.

As we started from an existing solver that we developped already for one NLS equation, to solve the system (1), we developped an algorithm that substitutes the solution of the first equation into the second equation. The numerical approximation is based on Crank-Nicholson scheme for the temporal discretization and Lagrange P1-Galerkin finite element method for spatial discretization.

By conformal-finite element approximation in H 1 0 (I), for I =] -l, l[ an open interval, N a positive integer, we then consider h = 1/(N + 1) a mesh space and the following approximation space:

V h = {v h ∈ C 0 ( Ī), v h (±l) = 0, v h|]-l+(i-1)h,-l+ih[ ∈ C[X], i = 1, . . . , N + 1}.
Obviously, the approximation space is spanned by classical hat functions denoted by (Ψ i ) 1≤i≤N and the nonlinearity is approached in a common way by:

-l+ih -l+(i-1)h |u| 2 uΨ k dx ≈ N j=1 |u j | 2 u j -l+ih -l+(i-1)h Ψ j (x)Ψ k (x)dx.
Its treatment has been done by a fixed-point iterative strategy which was not so coslty. Let us consider δt the time step, and mention by an upper index n the approximate value of the component of the solution at time t n = n δt, then

u(t n , .) ∈ V h , that means that ∃u n j , u(t n , x) = N j=1 u n j Ψ j (x), x ∈ I.
To solve the first equation, the numerical procedure consists of giving an initial guest u n+1,0 j as u n j and solving the following systems as soon as the threshold is not reached, u n+1,l+1 j -u n+1,l j > ε:

N j=1 u n+1,l+1 j Ψ j , Ψ k + ıδt 2 N j=1 u n+1,l+1 j ∇Ψ j , ∇Ψ k - γδt 2 N j=1 u n+1,l+1 j Ψ j , Ψ k -ıδt N j=1 (g 11 |u n+1,l j | 2 + g 12 |v n j | 2 )u n+1,l+1 j Ψ j , Ψ k = γδt 2 N j=1 u n j Ψ j , Ψ k - ıδt 2 N j=1 u n j ∇Ψ j , ∇Ψ k + N j=1 u n j Ψ j , Ψ k - ıκδt 2 N j=1 v n j ∇Ψ j , ∇Ψ k .
At convergence, one obtains u n+1 j 1≤j≤N whose values are substituting into the second equation.

A computer program has been developped in Python interfacing with Fortran that solves the generalized Manakov system (1) in a finite domain ] -l, l[ endowed by homogeneous Dirichlet boundary continions, with initial conditions (at t = 0) that have been given by Gaussian beams:

u 0 (x) = A π 1/4 a 1/2 exp - x 2 2a 2 , v 0 (x) = B π 1/4 b 1/2 exp - x 2 2b 2 .
In a first test series, we report the results for the set of initial conditions and parameters that have considered by other authors with some more details ( [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF]).

The computational domain has been taken sufficiently large that the numerical value of the initial conditions nearly verifies the homogeneous Dirichlet boundary conditions at the boundary locations. The first tests serie stands for different initial conditions (different values of (A, a, B, b)) with different γ and κ (such that γ < κ) for the case of the Manakov nonlinearity i.e. g 11 = g 22 = g 12 = 1 (see Fig. 1). In the case of Manakov nonlinearity, the results are in agreement of our theoretical computations. The third pattern of Fig. 1 indicates that the density oscillates in time and that the semi-norm of the solution stays bounded. Also, from the second pattern, one recovers numerically that the norm of each component of the solution oscillates but not in phase.

Let us recall that one property already known, is that

S 1 = I R (ūv + uv)dx,
the first Stokes integral of the solution of the Cauchy problem related to the Manakov system must stay constant in time [START_REF] Pelinovsky | Global Existence of solutions to coupled PT -symmetric nonlinear Schrödinger equations[END_REF]. The third pattern of Fig. 1 also numerically proves that the solution that is obtained has a Stokes integral S 1 which stays constant in time. The second tests serie also tends to illustrate that the H 1 -norm of the solution of the Cauchy problem for the PT-symmetric system (1) can blow up in finite time.

Conclusions

In this work we confirm the result that has been found previously by other authors. Despite the fact that we were not able to clarify if the H 1 -norm of the solution of the generalized Manakov system is globally bounded, we improved the estimate of the solution. Numerical results are also given to illustrate the theoretical result with a suitable numerical procedure.

Figure 1 :

 1 Figure 1: Numerical results of the Cauchy problem for the PT-symmetric system (1) in Manakov case (g 11 = g 22 = g 12 = 1). For each line the first plot stands for the initial condition, the second one traces the L 2 -norm of u, v and its gradient. The last graph represents the behavior of the density Q, the semi-norm D of the solution and one Stokes integral S 1 .

Figure 2 :

 2 Figure 2: Numerical results of the Cauchy problem for the PT-symmetric system (1) in non-Manakov case (g 11 = g 22 = 1 and g 12 = 1/2). For each line the first plot stands for the initial condition, the second one traces the L 2 -norm of u, v and its gradient. The last graph represents the behavior of the density Q, the semi-norm D of the solution and one Stokes integral S 1 .

  2 L 2 and v 2 L 2 according to the values of parameters. But our concern is to derive an upper bound which is valid for each case, if the sign of u 2 L 2 t and v 2 L 2 t are the same or not. Let us give details in the case of the sign of each component is different and that their sign change from an interval of the subdivision to the next one. For instance, (|u| 2 ) t ≤ 0 and (|v| 2 ) t ≥ 0 in [t 1 , t 2 ] whereas (|u| 2 ) t ≥ 0 and (|v| 2 ) t ≤ 0 in [t 2 , t 3 ]. Let us denote by I, the integral that follows:

	0	T	I R	arctan	Re u Im u	|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt
		≤	t 1 t 0	I R		arctan	Re u Im u	|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt+
	N -1 j=0	t 2 t 1	I R		arctan	Re u Im u		|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt
		+		t 3 t 2	I R		arctan	Re u Im u		|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt .
	Thus, by identyfying the sign of each term, and denoting the next left hand
	side by I j
	t 2 t 1	I R	arctan	Re u Im u	|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt
		+	t 3 t 2	I R	arctan		Re u Im u		|u| 2	t	+ arctan	Re v Im v	|v| 2	t	dx dt
	≤	t 2 t 1	I R	arctan	Re v Im v		|v| 2	t	-arctan	Re u Im u	|u| 2	t	dx dt
		+	t 3 t 2	I R		arctan		Re u Im u		|u| 2	t	-arctan	Re v Im v	|v| 2	t	dx dt

Acknowledgement

The authors wish to thank J. Laminie and E. Zahrouni for helpful comments and discussions.


Alina Quereilhac

Thierry Turletti

Walid Dabbous

Managing heterogeneous ns-3 experiments with NEPI

teaching and research institutions in France or abroad, or from public or private research centers.

 Objectives of NEPI

 Support all stages in the experiment life-cycle.  Automate repetitive and complex tasks.  Provide a platform-independent API.  Make it easy to conduct experiments that mix heterogeneous platforms (simulators, emulators, physical testbeds).

I -NEPI



  What is NEPI? The Network Experimentation Programming Interface (NEPI), is an experiment life-cycle management tool that enables conducting heterogeneous experiments where PlanetLab nodes and Emulated nodes can be seamlessly integrated with ns-3 simulations.



  The Network Experimentation Frontend (NEF), is a Graphical User Interface for NEPI that allows to design experiments by drag & dropping and interconnecting boxes in a canvas.  Deployment is done simply by selecting an option on the main Menu.  NEF provides several runtime views to monitor the evolution of the experiment, such as: o Topology View, to change experiment configuration during runtime. o Mobility View, to view the 2D changes in mobile nodes o Results View, to list and retrieve experiment results IV -Graphical User Interface  NEPI allows to build scenarios where ns-3 simulated networks can be integrated with the PlanetLab network. This is achieved by connecting an edge node inside the simulation to a virtual interface in the PlanetLab node where the simulation is running. The traffic generated inside the simulation can then be tunneled across the Internet to reach other PlanetLab nodes.  ns-3 simulations can also be integrated with emulated Ethernet networks, making it easy to inject in the simulation traffic generated by arbitrary unmodified applications. III -Heterogeneous experiments  Stages in the life-cycle  Experiment design  An experiment description (ED) is created by connecting boxes, representing experiment components, to other boxes on named ports called connectors. Each connector defines specific rules on how to create, configure and connect components during experiment deployment.  Experiment execution  The ED generated during design is used as the input for the Experiment Controller (EC). The EC is responsible for supervising the whole experiment, and instantiating the Testbed Controllers (TC) in charge of supervising the different experimentation platform instances involved in the experiment.  Upon creation, the EC translates the ED into platform-independent messages (e.g. "create Node N"), and sends them to the corresponding TCs. When a TC receives a message, it performs a platform-specific action in response (e.g. a ns-3 TC would create a ns3::Node C++ object in the ns-3 simulation upon receiving a create node message).II -Life-cycle supportPlanète Project-Team: http://planete.inria.fr/ ns-3 project: http://www.nsnam.org/