Benjamin Elie
email: benjamin.elie@inria.fr

Yves Laprie

Simulating alveolar trills using a two-mass model of the tongue tip

Keywords: numbers: 43, 70, Bk, 43, 72, Ja

HAL is a

I. INTRODUCTION

During the production of trill consonants (bilabial, alveolar, or uvular), aerodynamic forces applied to a particular articulator (lips, tongue tip, and uvula, respectively) are known to be responsible for the oscillation of the target articulator. Its oscillation leads to a quasiperiodic amplitude modulation of the acoustic pressure waveform radiated at the lips. This quasi-periodic amplitude modulation is the main characteristic feature of trills. This paper focuses on the mechanical and acoustic modeling of the production of alveolar trills, which are produced by the oscillation of the tongue tip in the alveolar region. This trill is the most common across languages [START_REF] Maddieson | Patterns of sounds[END_REF] , in comparison with other places of articulation. However, this is still a hard-to-produce sound that is mastered late in the acquisition process [START_REF] Jimenez | Acquisition of Spanish consonants in children aged 3-5 years, 7 months[END_REF] . This suggests that it requires specific articulatory and aerodynamic conditions, which have been reported to be the reason of assimilation in neighboring phonetic segments [START_REF] Solé | Production requirements of apical trills and assimilatory behavior[END_REF] . Curiously, only a few studies have investigated the acoustic and aerodynamic characteristics of the production of alveolar trills. Acoustic properties have been deeply detailed via analysis of the speech signal [START_REF] Mendoza | Acoustic characteristics of trill productions by groups of Spanish children[END_REF] , the intraoral pressure 5 , and the effect on the glottal source 6 . Linguopalatal contacts, along with the speech signal characteristics have also been investigated [START_REF] Recasens | A study of /R/ and /r/ in the light of the "DAC" coarticulation model[END_REF] . Despite these works, there is still a lack of acoustic and/or mechanical models of the production of trills. To the best of our knowledge, the work by McGowan 8 is the sole attempt to study the production of trills using numerical simulations. He investigated the mechanisms involved in the self-oscillations of the tongue tip with a very simplified vocal tract. Simulations show that the oscillations of the tongue tip are sustained thanks to the wall compliance of the vocal tract, hence a behavior similar to a one-mass system. The author also briefly discusses the potential impact of classic two-mass lumped models, similar to those used for the vocal folds [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF][START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] . Indeed, two-mass models (2MM) present the advantage of simulating phase shifts between the upstream and downstream parts of the vocal folds, as observed in real world experiments [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF] , and also to take a mobile flow separation point [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF][START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified two-mass model[END_REF] into account. One might expect a similar behavior for the tongue tip. The reader may find a detailed review of existing mass-spring systems for modeling the self-oscillation of the vocal folds in Erath et al. [START_REF] Erath | A review of lumped-element models of voiced speech[END_REF] , for instance.

In this work, the glottal source is connected to the classic transmission line paradigm [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] , considering recent improvements [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] in order to deal with a modified self-oscillating model of the vocal folds to account for partial glottal abduction. This is one of the major contributions of the paper in comparison with the work of McGowan 8 , since it enables the influence of the acoustic pressure variations upstream of the lingual constriction due to the vocal folds motion to be investigated with real anatomic data. Area functions have been derived from real-time cineMRI using a sparse reconstruction technique [START_REF] Elie | High spatiotemporal cineMRI films using compressed sensing for acquiring articulatory data[END_REF] , detailed in Sec. II.

The other contribution is the possibility of the presented simulation framework to consider the incomplete closure of the vocal tract during linguopalatal contacts. Indeed, results from Solé 5 suggest that there is no complete occlusion during the production of consonant trills, since the low frequency component of the mouth airflow is constantly strictly positive. Consequently, when dealing with acoustic simulations of the alveolar trills, one should guarantee that tongue tip oscillations with, or without contacts can be taken into account, and also that linguopalatal contacts do not completely stop the main air path in the vocal tract. In Sec. IV C, the paper presents a solution to deal with this situation by adding a lateral acoustic waveguide connected both at the upstream and the downstream parts of the tongue. This lateral connection is made possible thanks to the recent Extended Single-Matrix Formulation (ESMF) of the vocal tract [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] . The simulation framework including our self-oscillating tongue tip model is finally used to investigate the influence of various parameters on the production of alveolar trills. Results of the simulations are shown and discussed in Sec. V.

II. ACOUSTIC AND ARTICULATORY CHARACTERISTICS OF ALVEOLAR TRILLS

A. Acoustic characteristics

Alveolar trills [r] are sounds that are characterized by a quasi-periodic oscillation of the tongue tip. The vibration of the articulator is not directly controlled by the speaker, namely it is not due to an active muscular control, but is due to the aerodynamic forces applied at their vicinity [START_REF] Mcgowan | Tongue-tip trills and vocal-tract wall compliance[END_REF] . Studies about the production of alveolar trills include the mastering process 2 by children, occurrence in world's languages 1 , coarticulatory effects [START_REF] Solé | Production requirements of apical trills and assimilatory behavior[END_REF][START_REF] Lewis | Coarticulatory effects on Spanish trill production[END_REF][START_REF] Recasens | Coarticulation in Catalan dark [l] and the alveolar trill: general implications for sound change[END_REF] , acoustic characteristics [START_REF] Mendoza | Acoustic characteristics of trill productions by groups of Spanish children[END_REF]5 , and air-tissue interaction modeling 8 .

From these studies, it appears that alveolar trills are difficult to produce 2 , and may interfer with neighboring phonetic segments 3 . Aerodynamic measurements on speakers show that a high intraoral pressure is required to produce trills, and even higher for voiceless alveolar trills 5 . The intraoral pressure show significant variability across subjects 5 . The oscillation frequency is in the range 28-33 Hz for voiceless trills and 26-29 Hz for voiced trills 5 . The open phase is 1.94 times longer than the closed phase for voiceless trills, and 1.34 times longer for voiced trills.

Electropalatographic data (EPG) show that linguopalatal contacts occur during the production of consonant trills [START_REF] Solé | Production requirements of apical trills and assimilatory behavior[END_REF][START_REF] Recasens | A study of /R/ and /r/ in the light of the "DAC" coarticulation model[END_REF][START_REF] Recasens | Coarticulation in Catalan dark [l] and the alveolar trill: general implications for sound change[END_REF] . However, contact patterns significantly vary across speakers, and even across repetitions of a single speaker. For a few speakers, contacts occur only at the edge of the tongue dorsum, and the alveolar region does not show any contact between the tongue tip and the hard palate. Hence, there is no evidence of systematic contacts between the tongue and the hard palate, or the teeth, during the oscillation cycle. It is not unlikely that the tongue tip freely vibrates slightly away from the alveolar zone. Additionally, in all studies [START_REF] Mendoza | Acoustic characteristics of trill productions by groups of Spanish children[END_REF]5 , the presented waveforms of the airflow measured at the mouth is not null during the so-called closed phase. This suggests that trills do not involve complete occlusion. It is also clear in Fig. 3 of the paper of Solé 5 that the mouth airflow is systematically strictly positive during the trill production.

Consequently, when dealing with numerical simulations, one should consider the possibility of taking into account the incomplete closure of the vocal tract during linguopalatal contacts.

Fig. 1 shows two examples of acoustic speech signal of alveolar trills uttered in the pseudoword /ara/. The time envelope of the acoustic pressure waveform shows quasi-periodic cycles during the production of the alveolar trill [r]. It oscillates between a maximum, at a level slightly less than the adjacent vowels, and a minimum, which is not null.

B. Articulatory characteristics

Observing tongue tip movements during alveolar trills is very difficult. Indeed, since the oscillation has a small amplitude and is relatively fast, it requires high spatiotemporal resolutions. Thus, high-speed Ultrasound imaging and Electromagnetic articulography (EMA) 18 have been used to observe the time evolution of the tongue tip position during alveolar trills.

These studies provide data from alveolar trill realizations of German 18 , Czech 19 , Spanish, The simulation framework presented in this study requires the whole area function of the vocal tract. For that purpose, the vocal tract geometry has been observed via real-time cineMRI techniques. It consists in a reconstruction of the midsagittal view of the vocal tract at a reasonably high frame rate, using sparse reconstruction methods [START_REF] Elie | High spatiotemporal cineMRI films using compressed sensing for acquiring articulatory data[END_REF] . The frame rate is 29 frames per second, with a spatial resolution of 1×1 mm 2 . MRI experiments were performed on a 3T Signa HDxt MR system (General Electric Healthcare, Milwaukee, WI). Speech MRI data were obtained from a healthy volunteer with written informed consent and approval of local ethics committee. He is a 35 years old male and a French native speaker.

Although French does not have alveolar trills, our speaker's ability of uttering alveolar trills comes from a second language acquisition (Spanish) at an intermediate fluency level. The data were collected with a 16-channel neurovascular coil array. The protocol consisted in a sagittal slice through the middle of the vocal tract acquired with a custom modified Spoiled Fast Gradient Echo (Fast SPGR, TR 3.5ms, TE 1.1ms, line BW 83.33 kHz, flip angle 30 degrees, half echo in frequency direction, matrix 256×256, 512 temporal frames). For each temporal frame the modification consisted in acquiring only a randomized fixed size sample (10 lines) of the 256 phase lines of the k-space. The alveolar trill is pronounced in the intervocalic context /ara/.

We considered that acquiring MR images at 29 frames per second, namely similar to the trill frequency, is sufficient to get a good idea of the global position of the tongue tip during oscillations along with the rest of the vocal tract since we acquire several images per /r/ segment. For the study, the speaker was asked to slightly lengthen the duration of the /r/ segment, so that between 5 and 10 images per phonetic segment are available. Fig. 2 shows an example of the midsagittal view of the vocal tract during the production of an alveolar trill in the pseudoword /ara/. All the contours which are sufficiently invariant because they correspond to a non-deformable structure are tracked by correlation with a reference image chosen in the first images of the sequence. This is the case of the incisor roots for locating the mandible, and the upper part of the head and the vertebra used for locating the palate and pharyngeal walls. Deformable contours presenting a small number of deformation modes, like lips and larynx are tracked by means of a semi automatic method inspired from Fontecave and Berthommier [START_REF] Jallon | A semi-automatic method for extracting vocal-tract movements from X-ray films[END_REF] . The main principle is to delineate the target contour in key images randomly selected in the cineMRI, and then to index images where the contour has to be tracked with respect to the key images by using a Discrete Cosine Transform (DCT) distance. The image used for indexing images is limited to the region where the contour to be tracked is located. The resulting contour is obtained by averaging contours of the closest images. This semi-automatic tracking requires supervision from the user who can add key images to correct errors corresponding to an image too far from the existing key images. Finally, the tongue contour which presents a high variability, due to its movements and the dynamic MRI acquisition technique as well, has been delineated by hand.

From the contour extraction, area functions are obtained by, first, dividing the vocal tract shape in tubelets perpendicular to the centerline, and then applying α β transformations to recover the area [START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF] . This transformation is a power function that gives the cross-sectional area a(x) as a function of the midsagittal distance d(x):

a(x) = α(x)d(x) β(x) , (1)
where α and β are ad hoc coefficients. Their values depend on the region of the vocal tract, and x is the distance to the glottis along the centerline. Several studies have provided different values for α and β. We kept those given by Soquet et al. [START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF] for male speakers. The determination of the centerline plays a critical role in synthesis since it influences the length of the vocal tract. A specific algorithm was designed purposely [START_REF] Laprie | Articulatory copy synthesis from cine X-ray films[END_REF] .

As previously explained, the real-time cineMRI reconstruction enables several area functions to be recovered for a single [r] realization. Since the main aim of the paper is to investigate the possibility to simulate self-oscillations of the tongue tip in alveolar trills with lumped mass-spring models integrated into realistic shapes of the vocal tract, using area functions of alveolar trills on different contexts is not in the scope of the paper. Consequently, for the sake of brevity, we chose to use only one area function for the rest of the paper: it corresponds to the median area function extracted from the acquired sequence of images at the trill instants. The simultaneous audio recordings of the utterances allows each reconstructed MR image to be labeled with its corresponding produced phoneme. The MRI noise is removed from the recorded audio signal thanks to a source separation technique [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] .

The images corresponding to the trill instants are then the ones that correspond to the /r/ segment. The choice of a median area function is mainly for the purpose of compensating for the potential errors due to the manual delineation, as well as the variations in the alveolar constriction region due to the tongue tip movement during the trill production. The resulting area function is shown in Fig. 2. In the simulations, it corresponds to that of the vocal tract when the tongue is at its equilibrium position. Oscillations of the tongue tip will then modify the cross-sectional area of the lingual constriction, namely the area of the tubelet having the smallest area (at 15.9 cm in Fig. 2 (c)). Note that the back of the tongue creates a narrowing at the middle of the vocal tract. This has been previously reported as a common features of rhotic sounds, including alveolar and uvular trills [START_REF] Boyce | Acquiring rhoticity across languages: An ultrasound study of differentiating tongue movements[END_REF] .

III. PHYSICAL MODELS FOR SIMULATING THE TRILL PRODUCTION

The simulation framework used to compute the acoustic propagation inside the vocal tract is derived from the transmission line circuit analog model (TLCA) approach [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] . It considers plane waves propagating along a spatially sampled vocal tract, modeled as a set of connected acoustic tubes, or tubelets. Unlike the other widely used approach, the reflection type line analog (RTLA) model [START_REF] Kelly | Speech synthesis[END_REF][START_REF] Story | Phrase-level speech simulation with an airway modulation model of speech production[END_REF] , it easily deals with time-varying lengths of the vocal tract, and also with uneven spatial sampling of the vocal tract. For further information, the reader may find a detailed review of existing techniques for speech synthesis in Ref. [START_REF] Kröger | Articulatory synthesis of speech and singing: State of the art and suggestions for future research[END_REF] . The framework that is used in this paper, called ESMF 14 , considers recent improvements of TLCA-based techniques, such as the possibility to connect self-oscillating models of the vocal folds with a membranous glottal chink (a gap in the membranous portion of the glottis due to its partial abduction). In this paper, the acoustic model is extended to include the potential oscillations of the tongue tip.

A. Acoustic propagation

The vocal tract is seen as a waveguide network, in which each waveguide represents a side cavity [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF][START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches, Speech Communication 50(3)[END_REF] . In our study, side cavities include the sublingual cavity (cf. Fig. 2) and the lateral channels around the sides of the tongue tip. It has been shown, following the recent ESMF paradigm [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] , that the wave propagation inside such networks, connected to a self-oscillating model of the vocal folds, is driven by the following system of equations

f = Zu Z + Qu Q , (2)
where

f ∈ R (N +1) is a vector containing pressure forces, Z ∈ R (N +1)×(N +1
) is a tridiagonal matrix containing impedance and loss terms associated to each tubelet, Q is a square matrix the same size as Z having only one non-zero element, that is

Q (1,1) = R b , u Z ∈ R N +1
is the vector containing the volume velocities inside each tubelet, and

u Q ∈ R (N +1) = [U 2 1 , U 2 2 , . . . , U 2
N +1 ] T is the vector containing the square power of the volume velocities. The term R b is the Bernoulli resistance, as defined in Ref. [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] . Eq. 2 is computed at each time step, and the volume velocities u are the solution of the equation. The acoustic pressure P Out is computed as the first time derivative of the volume velocities at the lip termination, i.e. U N +1 , where N is the number of tubelets that models the vocal tract.

The acoustic propagation model supports the connection of self-oscillating models of the vocal folds with a membranous glottal chink [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] . In such case, the membranous portion of the glottis is composed of two distinct parts along the length of the vocal folds: an adducted part, modeled with a classic two-mass self-oscillating system [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] , and an abducted portion, i.e. the so-called membranous glottal chink, modeled as an acoustic branch. This allows incomplete closure of the glottis during oscillation cycles of the vocal folds due to the partial abduction of the glottis to be modeled, and has been found to be important for simulating breathy voice 30 , and voiced fricatives [START_REF] Elie | A glottal chink model for the synthesis of voiced fricatives[END_REF] . In this paper, the degree of abduction, denoted by D ab , is defined as the ratio, in percent, between the area of the membranous chink and the area of the fully abducted glottis. Consequently, D ab = 0% when there is no membranous chink, and D ab = 100% when the glottis is fully abducted.

B. Self-oscillating model of the tongue tip

In order to reproduce self-oscillations of the tongue tip, our model considers the lower wall of the tubelet corresponding to the lingual constriction to be a lumped mass-spring system, namely the 2MM, derived from the classic two-mass mechanical system for the vocal folds by Ishizaka and Flanagan [START_REF] Ishizaka | Synthesis of voiced sounds from a two-mass model of the vocal cords[END_REF] . In this paper, the self-oscillation is used to update the crosssectional area of the tubelet that corresponds to the lingual constriction, namely the tubelet that presents the minimal cross-section area. Consequently, only the extremity of the tongue tip is assumed to oscillate, without modifying the cross-section areas of the neighboring tubelets. In the following presentation of the model, the upstream and the downstream sections, as in Fig. 3, represent the tubelets located just downstream and upstream of the tubelet with the minimal cross-section area.

For simulating the self-oscillation of the vocal folds, the choice of a two-mass system instead of a single-mass one is motivated by the possibility to reproduce higher modes of oscillation, and especially phase shifts between the upstream and the downstream parts of the vocal folds. Indeed, at least two modes of oscillation have been extensively observed and shown to be important to consider in the simulations [START_REF] Titze | The physics of small-amplitude oscillation of the vocal folds[END_REF] . For tongue tip oscillations, this has been previously discussed by McGowan 8 , but there is, so far, no evidence of such a phase shift between the upstream and the downstream parts of the tongue tip during trill productions.

The 2MM used for the tongue tip is then similar to the one used for vocal folds. It includes continuous contours of the glottal constriction profile, a quasi-static mobility criterion for the flow separation point [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF][START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified two-mass model[END_REF] , and the computed pressure forces consider viscous losses [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] .

The geometry of the tongue tip is represented in Fig. 3. In this paper, indices u and d denote the upstream and downstream parts respectively. The position along the x-axis inside the lingual constriction is denoted by x i , with x = 0, 1, 2, 3, corresponding to the constriction input, the position of the upstream mass, the position of the downstream mass, and the constriction output, respectively. The total length of the lingual constriction is then l t = x 3 -x 0 , and its width is denoted w t . The equations for the tongue motion are where F u and F d are the pressure forces applied at the upstream and at the downstream masses, respectively. They are are defined as

     m u ḧu (t) + r u ḣu (t) + (k u + k c )∆h u (t) -k c ∆h d (t) = F u (t) m d ḧd (t) + r d ḣd (t) + (k d + k c )∆h d (t) -k c ∆h u (t) = F d (t) , (3)
F u (t) = w t x 1 x 0 x -x 0 x 1 -x 0 P (x, t)dx + w t x 2 x 1 x -x 2 x 1 -x 2 P (x, t)dx F d (t) = w t x 2 x 1 x -x 1 x 2 -x 1 P (x, t)dx + w t x 3 x 2 x -x 3 x 2 -x 3 P (x, t)dx,
To compute the pressure distribution P c (x, t) inside the lingual constriction, the airflow in the oral channel is considered as incompressible. Additionally, the viscous effects, which cannot be neglected when the lingual constriction is almost closed, are taken into account via a Poiseuille corrective term. Our model considers the viscous effects as predominant over the inertial effects so that the latter are not taken into account, as suggested by Lous et al [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] and Deverge et al. [START_REF] Deverge | Influence of collision on the flow through in-vitro rigid models of the vocal folds[END_REF] . Considering these assumptions, the pressure distribution P c (x, t) is given by

P c (x, t) = P u (x, t) + Be(x, t) + P o(x, t) (4)
where Be(x, t) and P o(x, t) are respectively the steady term of the Bernoulli equation and the Poiseuille corrective term. They are defined as:

Be(x, t) = - ρU 2 t (t) 2w 2 t 1 h 2 t (x, t) - 1 h 2 t (x 0 , t) , P o(x, t) = - 12µU t (t) w t x x 0 dx h 3 t (x, t) , (5)
where U t (t) is the instantaneous volume velocity of the airflow inside the lingual constriction.

Eq. ( 5) is valid for x located upstream of the separation point. Indeed, in this paper, similarly to the vocal folds model, the two-mass model of the tongue tip includes a mobile separation point x s . For the vocal folds model, the separation point has been shown to have a significant impact on the oscillating behavior [START_REF] Sciamarella | On the acoustic sensitivity of a symmetrical two-mass model of the vocal folds to the variation of control parameters[END_REF][START_REF] Zhang | Influence of flow separation location on phonation onset[END_REF] . However, predicting its location accurately is a very hard task as it requires numerical simulations of the full three-dimensional flow [START_REF] Zhang | Influence of flow separation location on phonation onset[END_REF][START_REF] Zhao | Computational aeroacoustics of phonation, Part I: Computational methods and sound generation mechanisms[END_REF] .

For lumped models as used in this paper, simplified methods to predict the location of the separation point have been proposed [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF][START_REF] Pelorson | Theoretical and experimental study of quasisteady-flow separation within the glottis during phonation. Application to a modified two-mass model[END_REF] . They usually consider the separation point to be located at a point downstream of the minimal glottal constriction, where the cross-section area is that of the minimal glottal constriction multiplied by an ad hoc criterion, usually taken between 1.15, as in Haas et al. [START_REF] Haas | Study of the effect of a moderate asymmetry on a replica of the vocal folds[END_REF] and 1.2, as in Lous et al. [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] . This is a very strong assumption, which may lead to inaccurate predictions of the separation point location, as it has been shown by numerical studies [START_REF] Sciamarella | Solving for unsteady airflow in a glottal model with immersed moving boundaries[END_REF] . Yet, this paper considers the quasi-static criterion as there is still a lack of accurate prediction methods suitable to lumped models.

Thus, considering the ad hoc criterion proposed by Lous et al. [START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] of 1.2, the separation The parameters chosen for the 2MM are summarized in Tab. I. The width of the lingual constriction w t and the height at rest h 0 are derived from the MRI data. Since lumped massspring systems models consider the shape of the constriction as rectangular, w t is estimated such that the constriction area a t = h 0 w t = α t h β 0 matches the cross-sectional area given by the power transformation of Eq. ( 1), hence

point x s is defined such that h s = h(x s ) = 1.2h u if 1.2h u < h d ,
w t = α t h βt-1 0 , (6)
where α t and β t are the transformation coefficients from Soquet et al. [START_REF] Soquet | Mid-sagittal cut to area function transformations: Direct measurements of mid-sagittal distance and area with MRI[END_REF] in the alveolar region, namely α t = 1.92 and β t = 1.2 when w t and h 0 are expressed in cm. The stiffness Production of alveolar trills (Version produced by the authors, accepted for publication in JASA)

values have been chosen according to the mass values so that the resonance frequency of the uncoupled mechanical system is 25 Hz, which is in the range of typical values observed in real speakers 5 . Assuming the mass density of the tongue ρ t is slightly more than water 39 , i.e. ρ t 1.04×10 6 g.m -3 , the distance between the base and the extremity of the tongue tip h tt , measured in MR images, is around 1.1 cm, and the tongue tip is a parallelepiped, each mass is

m = 1 2 ρ t × l t × h tt × w t 0.26g.
The total mass of the system is then 0.

B. Simulation of an alveolar trill with the lingual constriction modeled as a single waveguide

In this section, the lingual constriction is seen as a single waveguide. Results of the simulations are displayed in Fig. 4. It shows the oscillation of the lingual constriction height around its equilibrium position. In this example, quasi-periodic contacts between the tongue tip and the hard palate occur. The time envelope of the acoustic pressure is then a quasiperiodic function, which is at its minimum when the tongue tip is close to the hard palate (or is in contact with), and at its maximum when the tongue tip is far from the hard palate.

The displacement of the upstream and downstream masses of the tongue tip model shows no significant phase shift, which suggests a behavior similar to a single-mass system. This is an important result to consider for future studies since it implies that our flow model, which uses the quasi-steady flow approximation may be inadequate in the closing phase.

Indeed, Deverge et al. [START_REF] Deverge | Influence of collision on the flow through in-vitro rigid models of the vocal folds[END_REF] have shown that the flow unsteadiness is significant in the case of straight uniform channels. It is also worth noting that, in the case of almost straight uniform channels, the quasi-static mobility criterion for the flow separation point has no significant effect since the channel divergence is never larger than the chosen threshold. In order to address these important questions, improvements of the flow model should be considered in the future studies.

As reported by Dhananjaya et al. 6 , trills are known to have a significant impact on the glottal source. Our simulations also reveal that the tongue tip oscillations may disturb the motion of the vocal folds. Indeed, the amplitude of the oscillation of the vocal folds increases when the tongue tip is close to, or in contact with, the hard palate. This is visible in the bottom plot in Fig. 4, representing the glottal opening h g , where the maximal values of h g slightly increase during two glottal pulses at t = 8ms, t = 48ms, and t = 88 ms. For each of these time instants, the tongue tip is in contact with the hard palate. Basically, the amplitude of the glottal opening waveform follows the lingual constriction closure with a small time delay of a few milliseconds. More interestingly, the downstream part of the vocal folds is more impacted, as it exhibits a larger oscillation amplitude than the upstream part in the closing phase of the lingual constriction. Since the closing phase of alveolar trills is known to exhibit a rise of the supraglottal pressure 5 , this is certainly the main reason of this modification of the glottal opening waveform. The upstream part of the vocal folds is the most impacted since it is directly connected to supraglottal part of the vocal system.

However, the acoustic and the airflow waveforms P Out and U m show discrepancies with real alveolar trills. Indeed, the simulated waveform exhibits null values during the closed phase, while this is never observed in real trills. This is due to the fact that the vocal tract is seen here as a single waveguide, which leads to the complete occlusion of the vocal tract during linguopalatal contacts. The next section presents a method to fix this issue by adding lateral channels that allows incomplete closures to be simulated.

C. Simulations of the incomplete closure of the vocal tract at the linguopalatal contact instants

As discussed in Sec. II, in real trills linguopalatal contacts do not lead to the complete occlusion of the oral channel. Indeed, both the air flow at the mouth and the acoustic pressure waveform are not null. In order to simulate the incomplete closure of the vocal tract during linguopalatal contacts, we suggest to add a lateral acoustic branch to simulate the air path going along the side of the tongue. The ESMF paradigm that is used in this paper enables such situations [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] to be taken into account in the acoustic simulations. Due to the lack of data (cineMRI films provide only midsagittal views of the vocal tract), the area function of the lateral branch is estimated according to these assumptions: i) connections at the upstream and at the downstream parts are located respectively at 14.0 cm and 16.2 cm from the glottis, ii) the cross-section area of the lateral branch is taken at a certain ratio r l of the vocal tract area function located between the connections, iii) the length of the lateral branch is the same as the length of the main air path. The presented simulation uses the 2MM of the tongue tip with the same mechanical parameters as in Tab. I, integrated into a vocal tract modeled with a lateral branch that is connected to the main path with an area ratio r l of 50%. The simulation with a ratio r l = 50%, displayed in Fig. 5, shows that this vocal tract configuration enables linguopalatal contact with incomplete closure of the vocal tract to be simulated. When the tongue tip and the hard palate are in contact (h t ≤ 0), the air flow U m is not null, nor is the acoustic signal pressure. These features are in agreement with what is commonly observed in real trills (see Fig. 1, for instance).

Note that the inclusion of the lateral waveguide does not significantly modify the global behavior of the tongue tip oscillations, except for a slight decrease of the trill amplitude.

There is still no phase shift between the upstream and the downstream parts of the tongue tip, still suggesting a behavior similar to a single mass system. Note also that the rise of the vocal folds opening during closed phase is still visible at t = 10 ms, t = 50ms, and t = 90ms.

Production of alveolar trills (Version produced by the authors, accepted for publication in JASA)

V. CONDITIONS OF PRODUCTION OF ALVEOLAR TRILLS

This section investigates the influence of various configurations of the vocal tract and of the glottis on the simulation of the tongue tip oscillation. Studied parameters include the subglottal pressure, the mass of the tongue tip, the height of the lingual constriction at rest, the glottal abduction degree, and the lateralization of the air flow.

A. Trill characteristic features

The analyzed trill characteristic features are the trill frequency f t , the trill amplitude, denoted ĥt , and the trill contact ratio, denoted C r . For each simulation of sustained trills, a group of 10 successive periods are considered. Fig. 6 illustrates the identification of the different features. For each oscillating cycle n, A n is the difference between the maximal and the next minimal lingual constriction height, T n is the time difference between two successive maxima, and c n is a binary number, equal to 0 if there is no linguopalatal contact during cycle n, and 1 if there is a contact. The trill frequency is the inverse of the mean period value, namely Fig. 7 shows the trill contact ratio and the trill amplitude as a function of the values of each mass modeling the tongue tip and the subglottal pressure. Not surprisingly, the trill amplitude is the largest for small values of the tongue mass and for large values of the subglottal pressure P sub . For mass values smaller than 0.1 g, the amplitude of the tongue tip oscillation may be larger than 1 cm, which seems unrealistically high, since it is out of the range of reported data on real speakers 18,[START_REF] Howson | Examination of the grooving patterns of the Czech trill-fricative[END_REF][START_REF] Kim | Toward automatic vocal tract area function estimation from accelerated three-dimensional magnetic resonance imaging[END_REF] . This leads to area variations of approximately 1.4 cm 2 , according to Eq. ( 1), and is much larger than data provided by McGowan 8 , who reported maximal constriction area between 0.27 and 0.67 cm 2 . These areas correspond to maximal constriction heights of 0.44 cm and 0.21 cm, respectively. In the range of subglottal pressure values usually observed in natural speech [START_REF] Giovanni | Estimated subglottic pressure in normal and dysphonic subjects[END_REF] , i.e. under 1000 Pa, the trill amplitude values that lies in this range (0.21-0.44 cm) are reached for mass values between 0.2 g and 0.4 g. In this range, linguopalatal contacts occur almost systematically (C r > 80%). Note that the contact ratio C r follows the same evolution as the trill amplitude. For the rest of the paper, both masses modeling the tongue will be set to 0.26 g, as in Sec. IV, since it roughly corresponds to anatomical data, and additionally lies in the range of mass values that simulate realistic oscillations of the tongue tip.

f t = 10

C. Effect of the tongue initial position

For these simulations, the glottis is entirely adducted, i.e. D ab = 0%, and there is no lateralization, i.e. r l = 0%. The length and the width of the lingual constriction are set The tongue initial position h 0 is set to different values to study its impact on the selfoscillating motion of the tongue tip. Fig. 8 shows the trill frequency f t , the trill amplitude ĥt , and the trill contact ratio C r as a function of the subglottal pressure P sub and the initial tongue position h 0 . It clearly shows that the initial tongue position has a significant influence on the production of alveolar trills. For instance, for a given P sub , the trill amplitude as a function of h 0 exhibits a local maximum, i.e. when h 0 is between 0.5 and 1 mm. This implies that the tongue should not be in contact with the hard palate at the beginning of the trill.

Results also show that when h 0 is larger than approximately 1.25 mm, linguopalatal contacts occur (C r >= 10%) only when P sub is larger than a certain threshold, between 1400 and 1500 Pa, for h 0 ≤ 0.25cm, and larger than 1500 Pa for h 0 ≥ 0.25cm. This threshold value increases as h 0 increases: it is larger than 2000 Pa for h 0 ≥ 0.3cm. These observations show that the alveolar trills are favored when the tongue tip is very close to the hard palate, i.e. less than 1.25 mm, but not in contact. Consequently, producing alveolar trills requires a very fine adjustment of the tongue tip position in order to generate the self-oscillation; which confirms the fact that alveolar trills are difficult to produce. It is also worth noting that larger lingual constriction heights at rest tend to increase the trill frequency. The trill frequency is the most sensitive to the equilibrium position when h 0 is smaller than the local maximum, i.e. when h 0 < 0.5 mm.

D. Effect of the incomplete closure of the vocal tract

For these simulations, the glottis is entirely adducted, i.e. D ab = 0%, the equilibrium position, the length and the width of the lingual constriction are set to the values shown in Tab. I. The masses m 1 and m 2 are set to 0.26 g. The lateral area ratio r l varies from 0 to 50%, and the subglottal pressure from 500 to 2000 Pa.

The incomplete closure of the vocal tract is studied via the lateralization ratio r l . Fig. 9 shows that this has no significant effect on the trill characteristics: the trill frequency is barely lowered by 1 Hz when there is an incomplete closure of the vocal tract, and, for a given P sub , the trill amplitude and the trill contact ratio exhibit a slight tendency to decrease as the incomplete closure increases (larger r l ). However, although the incomplete closure of the vocal tract during linguopalatal contacts has a very small impact on the mechanical behavior of the tongue tip, it has been shown in Sec. IV C that the acoustic impact is not negligible. Consequently, when dealing with acoustic synthesis of alveolar trills using physical models, one should consider the lateralization of the airflow.

E. Effect of the glottal abduction

For these simulations, there is no lateralization, i.e. r l = 0%, the equilibrium position, the length and the width of the lingual constriction are set to the values shown in Tab. I.

The masses m 1 and m 2 are set to 0.26 g. The degree of glottal abduction D ab varies from 0 to 100%, and the subglottal pressure from 500 to 2000 Pa. Fig. 10 shows that the presence of the partial abduction has a significant impact on the trill characteristics. Indeed, it has for main effect, for a given P sub , to increase the trill frequency, the trill amplitude, and also the trill contact ratio. Indeed, when the degree of abduction D ab is greater than approximately 30%, the trill amplitude is smaller than 2 mm, and the trill contact ratio shows that there is almost no linguopalatal contacts.

There is a transition zone, when D ab is between 20 and 40%, where the trill contact ratio is very sensitive to small perturbations of the glottal abduction. Before this transition, for D ab < 25%, the trill contact ratio is 100%, and it vanishes for D ab > 40%, except when P sub is larger than 1700 pa, approximately.

Note that the increase of the trill frequency for voiceless trills is in agreement with the experimental observations by Solé 5 , who reported higher trill frequencies for voiceless trills (mean 29.3 Hz) than for voiced trills (mean 28.1 Hz). Another explanation of a higher trill frequency of voiceless trills could be the fact that they require higher subglottal pressures. Indeed, for a given abduction degree, increasing the subglottal pressure increases the trill frequency, as shown in Fig. 10 (a). The two-mass model of the tongue tip oscillates at frequencies similar to those observed in natural speech 5,8 , i.e. in the range 25-30 Hz. In our simulations, there are no significant phase shifts between the two masses, which suggests that, unlike the vocal folds, its behavior is very similar to the single-mass model.

VI. CONCLUSION AND FURTHER WORKS

The simulation framework presented in this paper also deals with the commonly observed case of the incomplete closure of the vocal tract during linguopalatal contacts. It is taken into account via the connection of lateral branches around the tongue. Although our simulations show that the connection of lateral branches does not significantly modify the oscillation characteristics of the tongue tip, it has a strong acoustic impact: the intraoral airflow and the acoustic pressure waveforms are not completely null in the closed phase of the trill, as it has been observed in natural speech in all of the previous studies 5, [START_REF] Mcgowan | Tongue-tip trills and vocal-tract wall compliance[END_REF] . Consequently, the incomplete closure of the vocal tract should be taken into account when dealing with the acoustic synthesis of alveolar trills using physical models.

Additionally, our simulations have revealed the fact that there is an initial position of the tongue tip that favors the production of alveolar trills. Indeed, for a given subglottal pressure, the trill amplitude as a function of the equilibrium position of the tongue tip admits a local maximum, located between 0.5 and 1 mm off the hard palate. When this tonguepalate distance at rest is greater than 1 mm, the trill amplitude significantly decreases, making the production of alveolar trills more difficult. This might be one of the reasons why consonant trills are hard-to-produce sounds in human spoken languages.

Finally, it has been shown that the amplitude of the tongue oscillation is significantly reduced in voiceless configurations, i.e. when the vocal folds are partially, or completely abducted. This is especially true when the degree of glottal abduction exceeds a critical value, between 30 and 40% of the total glottal abduction, depending on the value of the subglottal pressure. This suggests that voiceless trills require higher subglottal pressure, and/or additional adjustments, to produce trills with a similar oscillation amplitude than voiced trills. This is a possible explanation for the fact that voiceless trills are rather rare among human languages in comparison with voiced trills as they require more specific efforts.

This first attempt to model the tongue tip oscillations with classic lumped mass-spring systems opens new ways of investigations of the production of trills. Indeed, one could imagine similar models applied to other trills, such at the lips for the bilabial trills, or at the uvula for the uvular trills. The possibility to support realistic geometries of the vocal tract derived form cineMRI acquisitions of real speakers is a useful tool to study the impact of various articulatory configurations on the production of trills, and consequently to investigate the reasons of the difficulty encountered by numerous speakers to produce trills. Additional data, considering phonological contexts, are planned to be acquired in the near future to complete this investigation. Finally, the paper raises several questions about the dynamics of the tongue tip oscillation and about the aerodynamic behavior of the flow in the lingual constriction during alveolar trills that remain unsolved, mainly due to the lack of full numerical simulations, and/or fine experimental observations. Thus, this paper may constitute a base for such investigations in order to eventually address these important questions for the full understanding of the alveolar trill dynamics. More specifically, the questions about the flow separation point and the phase shift should be addressed with more realistic flow models, and with rounded geometries of the tongue tip that better approximate its real shape.

Figure 1 .

 1 Figure 1. Audio speech signal of two utterances of the pseudoword /ara/, with a zoomed-in version of the alveolar trill [r]. The time envelope is plotted as a dashed line

Figure 2 .

 2 Figure 2. (a) MR image of a midsagittal view of the vocal tract during the production of an alveolar trill in the sequence /ara/. Contours of the articulators are shown in the MR image, as well as (b) the spatial sampling of the vocal tract into acoustic tubelets. (c) resulting median area function obtained from real-time cineMRI data.

Figure 3 .

 3 Figure 3. Two-mass model for the tongue tip. P d , P u , and P c denote the pressure downstream, upstream and inside the lingual constriction, respectively. The mass, stiffness and the damping of the tongue are denoted by the triplet {m, k, r}. They are coupled by a the spring k c . The height of the lingual constriction is denoted by h.

. 7 Length of the lingual constriction l t mm 4. 4 Width of the lingual constriction w t mm 10 . 5 Mass m = m 1 = m 2 g 0. 26 Stiffness k = k 1 = k 2

 74105262 and x s = x d otherwise.Finally, P c (x, t) = P d (t) if x > x s .The collision model uses artificial stiffening of the tongue tip, namely the stiffness parameters k u and k d are multiplied by a factor κ if h u , or h d , are equal or less than 0, respectively, and the damping is multiplied by a factor η. These factors are also ad hoc criteria, and their estimation is discussed in Sec. IV B.IV. SIMULATIONS OF ALVEOLAR TRILLSUsing the area function derived from cineMRI data of a real speaker, sustained alveolar trills are simulated with the 2MM. First, an alveolar trill is simulated without a lateral waveguide around the lingual constriction, such that linguopalatal contacts create a complete Production of alveolar trills (Version produced by the authors, accepted for publication in JASA) occlusion of the vocal tract. Then, a solution consisting in adding a lateral branch around the lingual constriction is presented to simulate the incomplete closure of the vocal tract during linguopalatal contacts.A. Parameters of the tongue tip modelTable I. Parameters of the two-mass model of the tongue tip Parameter Unit Value Subglottal pressure P sub Pa 1000 Equilibrium position of the tongue tip h 0 mm 0

  26 × 2 = 0.52g. The ad hoc criteria κ and η of the collision model are estimated empirically. Their values impact the closed phase duration, and consequently, the open quotient of the tongue tip oscillation, namely the ratio of the duration of the open phase in comparison with the trill period. They have been estimated as the values of {κ, η} such that the framework simulates trills with an open quotient similar to those experimentally estimated on real speakers, namely 57.3 %, as reported by Solé 5 . This open quotient has been reached for κ = 1.5 and η = 1.2. The collision model uses these values for the rest of the paper.

Figure 4 .

 4 Figure 4. Simulation of an alveolar trill with linguopalatal contacts using the two-mass model. From top to bottom are plotted the acoustic pressure radiated at the lips (P Out ) and its time envelope (dot-dashed line), the time evolution of the lingual constriction height h t at the upstream (solid line) and downstream (dashed line) sections, the low frequency filtered airflow at the mouth U m . The bottom figure plots the displacement h g of the upstream (solid line) and downstream (dashed line) parts of the vocal folds. Positive values of h t and h g indicate open lingual and glottal constrictions, respectively, while null and negative values indicate closed constrictions.

Figure 5 .

 5 Figure 5. Results of the simulation of a sustained alveolar trill with a lateral configuration (r l = 50%). From top to bottom: acoustic pressure radiated at the lips P Out (solid line) and its time envelope (dash-dotted line), lingual constriction height h t at the upstream (solid line) and at the downstream (dashed line) parts, and air flow U m at the mouth. The bottom plot represents the displacement h g of the vocal folds at the upstream (solid line) and at the downstream (dashed line) sections. Positive values of h t and h g indicate open lingual and glottal constrictions, respectively, while null and negative values indicate closed constrictions.

Figure 6 .

 6 Figure 6. Example of a typical lingual constriction aperture waveform and identifying features A n , T n , and c n , for the computation of the trill characteristics.

10 n=1 T n .

 n The trill amplitude ĥt is the mean value of the difference between the maximal and the next minimal lingual constriction height in the group of 10 successive oscillatory cycles, henceThe trill contact ratio C r is the percentage of oscillatory cycles where there is a linguopalatal contact among 10 successive periods: of the mass and the subglottal pressure For these simulations, the glottis is entirely adducted, i.e. D ab = 0%, and there is no lateralization, i.e. r l = 0%. The other parameters, namely the equilibrium position h 0 , the length and the width of the lingual constriction, and the contact model parameters κ and η, are set to the values shown in Tab. I. The masses m 1 and m 2 vary from 10 mg to 1 g, and the subglottal pressure from 500 to 2000 Pa. For each simulation, m 1 = m 2 .

Figure 7 .

 7 Figure 7. (a) Trill amplitude and (b) trill contact ratio, as a function of the mass m 1 = m 2 and the subglottal pressure P sub . Solid lines represent the contour lines.

  This paper has discussed the possibility to model the self-oscillation of the tongue tip in alveolar trills via classic lumped mass-spring systems. Based on a model similar to the one used for vocal folds, the two-mass model[START_REF] Lous | A symetrical two-mass vocal-fold model coupled to vocal tract and trachea, with application to prothesis design[END_REF] , it has been shown that a quasi-periodic oscillation of the tongue tip could be reproduced to simulate the production of alveolar trills.In comparison with previous works 8 , the simulation framework presented in this paper can support physiologically realistic area functions of the vocal tract, and also the connection with a self-oscillating model of the vocal folds with a membranous glottal leakage to simulate voiced and voiceless alveolar trills. It can also account for the incomplete closure of the vocal tract during linguopalatal contacts due to the lateral leakage of the airflow around the tongue.

Figure 8 .

 8 Figure 8. (a) Trill frequency f t , (b) trill amplitude ĥt , and (c) trill contact ratio C r , as a function of the subglottal pressure P sub and the equilibrium position of the tongue h 0 . Solid lines represent the contour lines.

M.-J. Solé, Aerodynamic characteristics of trills and phonological patterning, J. of Phon.

N. Dhananjaya, B. Yegnanarayana, P. Bhaskararao, Acoustic analysis of trill sounds, J.

ACKNOWLEDGEMENTS

The work belongs to the project ArtSpeech, with the financial support of the French ANR (Agence Nationale de la Recherche). The authors are grateful to the laboratory IADI (Imagerie Adaptative Diagnostique et Interventionnelle) and the CIC-IT (Centre d'Intervention Clinique -Innovations Technologiques) of Nancy for the MRI acquisitions. articulography: Use of alternating magnetic fields for tracking movements of multiple points inside and outside the vocal tract, Brain and Language 31 (1) (1987) 26-35.


Nicolas Schnepf

Rémi Badonnel

Abdelkader Lahmadi

Stephan Merz

Automated Verification of Security Chains in Software-Defined Networks with Synaptic

Keywords: Security Management, Software-Defined Networking, Formal Verification

Software-defined networks provide new facilities for deploying security mechanisms dynamically. In particular, it is possible to build and adjust security chains to protect the infrastructures, by combining different security functions, such as firewalls, intrusion detection systems and services for preventing data leakage. It is important to ensure that these security chains, in view of their complexity and dynamics, are consistent and do not include security violations. We propose in this paper an automated strategy for supporting the verification of security chains in software-defined networks. It relies on an architecture integrating formal verification methods for checking both the control and data planes of these chains, before their deployment. We describe algorithms for translating specifications of security chains into formal models that can then be verified by Satisfiability Modulo Theories (SMT) solving or model checking. Our solution is prototyped as a package, named Synaptic, built as an extension of the Frenetic family of SDN programming languages. The performances of our approach are evaluated through extensive experimentations based on the CVC4, veriT, and nuXmv checkers.

I. INTRODUCTION

The growing dynamics and size of the Internet pose new challenges in terms of security management. An illustrative example can be given with the case of connected devices, such as smartphones and smart objects, even if our work is not restricted to this specific context. The multiplication of these devices is an important factor of Internet growth. It also contributes to a larger attack surface, such as the recent attack against the Dyn DNS service, which relied on a botnet of infected smart objects. Malicious applications targeting these devices are also increasing massively every year. Preventive security methods that consist in analyzing the applications before their publication on markets show their limitations. For instance last year, Kaspersky Labs estimated that 2,961,727 malwares were published on the official Google application market [START_REF] La Polla | A Survey on Security for Mobile Devices[END_REF], [START_REF] Faruki | Android Security, a Survey of Issues, Malware Penetrations and Defenses[END_REF]. Security mechanisms must be dynamically adjusted to these evolutive threats. The resources of devices in terms of memory, cpu, and battery are also often limited, which may make the local deployment of protective mechanisms more challenging.

In the meantime, the programmability offered by softwaredefined networks (SDN) provides new perspectives with respect to protection [START_REF] Feamster | The Road to SDN, an Intellectual History of Programmable Networks[END_REF], [START_REF] Feamster | Software-Defined Networks: Improving Network Management with SDN[END_REF]. First, security mechanisms such as firewalls, intrusion detection systems, and data leakage prevention services, can be virtualized in cloud infrastructures, and be offered as outsourced services. We will refer to them as security functions in the remainder of this paper. Second, these functions can then be combined to build elaborate security chains that are deployed and adjusted, in a dynamic manner, to protect devices and their associated applications. These chains can be fully outsourced in the network, or with a hybrid in-cloud/on-device deployment scheme. Softwaredefined networking separates the control plane from the data plane and interconnects them by an open interface, which may typically be the OpenFlow protocol. In our case, the data plane corresponds to the forwarding switches with the security functions that are currently deployed in the network, whereas the control plane describes the orchestration underlying automatic modifications of the configuration of the data plane. In this way, different security functions can be combined as a chain and deployed in the network. This dynamic chaining allows for an adaptive response to attacks. However, it is crucial to ensure that these security chains, in view of their complexity and dynamics, are consistent. Misconfigurations may reduce their efficiency, and introduce security violations.

We therefore study in this paper the use of automated techniques for verifying security chains in software-defined networks. Our approach relies on a software-defined architecture and integrates formal verification methods for checking both the control and data planes. We introduce algorithms for translating specifications of security chains into formal models that are then automatically verified by SMT solving and model checking. Our solution has been prototyped as a package, named Synaptic, built as an extension of the Frenetic family of SDN programming languages. The solution is evaluated through extensive experiments using the checkers CVC4 [START_REF] Barrett | CVC4[END_REF], veriT [START_REF] Bouton | veriT: An Open, Trustable and Efficient SMT-Solver[END_REF] and nuXmv [START_REF] Cavada | The nuXmv symbolic model checker[END_REF]. Our main contributions are: (i) the design of an automated verification strategy and its architecture, (ii) algorithms for translating specifications of security chains into formal models, (iii) the prototyping of the Synaptic package as an extension of the Frenetic language family, and (iv) the performance evaluation with different tools for SMT solving and model checking.

The remainder of this paper is organized as follows. Section II gives an overview of existing work in the area. Section III describes our automated verification strategy, its architecture and the supporting algorithms for generating formal models based on security chains. It then presents our prototype, named Synaptic, built as an extension package for the Frenetic language family. Section IV describes performance results obtained with different verification tools. Section V concludes and points out future research perspectives.

II. RELATED WORK

The increasing development of cloud infrastructures has spurred the virtualization and outsourcing of a large variety of network services, including for security purposes. Network programmability provides a support for combining them in order to elaborate service chains. For instance, Sherry et al. [START_REF] Sherry | Making Middleboxes Someone Else's Problem: Network Processing as a Cloud Service[END_REF] propose a solution for dynamically and transparently outsourcing middleboxes across several cloud providers using virtualization and different redirection mechanisms. Gibb et al. [START_REF] Gibb | Outsourcing Network Functionality[END_REF] present a similar approach where a cloud-based architecture is designed for outsourcing network functionalities using SDN. Regarding the chaining of such network functionalities, Qazi et al. introduces SIMPLE [START_REF] Qazi | SIMPLE-fying Middlebox Policy Enforcement using SDN[END_REF], a policy enforcement layer based on SDN and flow correlation for middlebox traffic steering. In the same vein, Fayazbakhsh et al. propose Flowtags [START_REF] Fayazbakhsh | Enforcing Network-Wide Policies in the Presence of Dynamic Middlebox Actions using FlowTags[END_REF], an architecture where middleboxes are extended to support OpenFlow and to deal with dynamic middlebox flow mangling. In the area of smart devices, Sapiol et al. [START_REF] Sapio | Per-user Policy Enforcement on Mobile Apps Through Network Functions Virtualization[END_REF] define a Network Function Virtualization (NFV) router to provide a per-user policy enforcement on mobile applications through service chaining. Previous work of our group [START_REF] Hurel | Towards Cloud Based Compositions of Security Functions for Mobile Devices[END_REF], [START_REF]Behavioral and Dynamic Security Functions Chaining for Android Devices[END_REF] proposed the design of security chains for protecting such devices, and showed their benefits. The dynamics and multiplication of elaborated chains require verification techniques to ensure their consistency.

There exists a substantial body of literature for supporting the verification of data planes in the area of softwaredefined networking. The purpose is typically to prevent misconfigurations that could lead to black holes and loops in the network. For instance, VeriCon [START_REF] Ball | Vericon: Towards Verifying Controller Programs in Software-Defined Networks[END_REF] is a language for specifying software-defined networking policies and is accompanied by a solution for checking whether a policy verifies invariants expressed using predicate logic. It does not address the specification of temporal aspects of the control logic. NICE [START_REF] Canini | A Nice Way to Test OpenFlow Applications[END_REF] is another example of SDN verification based on unit tests; however, tests cannot ensure complete coverage of the behavior of complex systems such as SDN policies. FlowChecker [START_REF] Al-Shaer | FlowChecker, Configuration Analysis and Verification of Federated OpenFlow Infrastructures[END_REF] represents the data plane as a binary decision diagram (BDD), whereas properties are expressed in computation tree logic (CTL). The complexity of the BDDbased modeling may constitute an obstacle to its large-scale deployment. VeriFlow [START_REF] Khurshid | Veri-Flow: Verifying Network-wide Invariants in Real Time[END_REF] also proposes support for the realtime verification of OpenFlow rules, but its low level of abstraction may make it cumbersome to verify large security chains. Finally, VeriSDN [START_REF] Kang | A Verification Method of SDN Firewall Applications[END_REF] is another example of control plane verification in the specific case of SDN-based firewalls, although the verification of more sophisticated security functions is considered. In a similar spirit, work by Shaer et al. [START_REF] Al-Shaer | Discovery of policy anomalies in distributed firewalls[END_REF] considers anomaly detection in distributed firewalls. Most of these approaches are focusing on the verification of the data plane, and may miss dynamic aspects of the chain operating on the control and data planes.

Decoupling the functional specification of security chains and their translation into low-level configuration rules is an important requirement for enabling their formal verification, which motivates the use of high-level network programming languages. The Frenetic language family for network programming [START_REF] Foster | Frenetic, a Network Programming Language[END_REF] counts among the most prominent approaches for verifying the control plane. It includes a language, called Pyretic, for specifying network configurations in Python; these configurations are then compiled into low-level rules [START_REF] Foster | Languages for Software-Defined Networks[END_REF]. Pyretic is complemented by an extension, called Kinetic, for describing control plane policies; Kinetic also offers formal verification techniques [START_REF] Kim | Kinetic: Verifiable Dynamic Network Control[END_REF]. Whereas this work is restricted to the verification of the control plane, we argue in favor of verifying both the security functions chaining of the data plane and the dynamics of the control plane in an integrated manner in order to prevent errors due to potential misconfigurations.

III. AUTOMATED VERIFICATION OF SECURITY CHAINS

We propose in this paper an automated strategy for verifying both the control and data planes of security chains in softwaredefined networks before their deployment or re-configuration. Our goal is to analyze both the consistency of the data plane, taking into account its dynamic evolutions specified by the control plane. The solution is based on a dedicated architecture, which includes a verifier, named Synaptic, capable of translating security chain specifications into formal models that can be verified by SMT solving or model checking.

A. Proposed architecture

The architecture is depicted in Figure 1, with the illustrative context of protecting smart devices, such as smartphones. However, our work is not limited to this specific use case, since security chains are used more widely. The figure highlights three main entities: (i) on the left, the smart device, with several running applications, to be protected with an integrated agent (which can typically be instantiated by an OpenFlow virtual switch for redirection purposes); (ii) in the middle, a cloud provider infrastructure hosting several security functions as well as a security manager, and (iii) on the right, the remote destinations interacting with the applications executed on the smart device. The security manager is supported by an SDN controller that orchestrates security chains, using the Pyretic network programming language. In addition to the Kinetic extension that provides functions for verifying the control plane (blue), it exploits our Synaptic checker for verifying the data plane (shown in red). Finally security functions can be deployed either on the device or in a cloud infrastructure (dashed lines).

When an application initiates a communication with a remote destination, all the messages from and to that application go through the agent (virtual switch) of the device. The switch may probe the SDN controller from the cloud provider in order to know how to redirect the related messages for security checks. Depending on the risks and the context, the security manager activates the appropriate security functions. By pushing the necessary SDN rules within the cloud provider network, the controller composes the security functions to build an appropriate chain and notifies the switch. The latter directs the incoming and outgoing traffic through the chain before forwarding it to the final destination. Security functions can either be hosted locally on the smart device or in the cloud. Compositions of security functions are designed by the security manager according to several factors, including the originating application, the remote destination, and the network properties. For instance, the security functions may include network or applicative firewalls, intrusion detection systems, and mechanisms for preventing data leakage. They are not limited to traffic analysis, but may also include functions analyzing the configuration of smart devices and their applications. The design of these chains and the choice of security functions is out of the scope of this paper: we focus on their automated verification. Our solution extends the Frenetic framework. Specifically, the Pyretic network programming language is used for describing data plane configurations. The language includes four basic rules to express a static filtering on the traffic, namely: identity to forward all packets without conditions, drop to block all packets without condition, match(x=y) to forward only packets for which field x has the value y, and modify(x=y) to forward all packets and set the field x to the value y. They can be composed using three different composition operators: a sequential operator ( ), a parallel operator (+) and a negation operator (∼). In addition, the Kinetic language, also part of the Frenetic family, is intended for specifying and verifying the control plane, as a finite state machine. More specifically, it makes use of the NuSMV model checker to verify whether the corresponding automaton satisfies properties written as formulas of the temporal logic CTL. However, it does not verify the data plane.

B. Extension of the Frenetic framework

We therefore define an extension of the Pyretic language capable of verifying the data plane. Our checker, that we call Synaptic, is based on a translation of the rules of this language into formal models. We consider two categories of formal models for supporting this verification. The first category corresponds to SMT solving. It is based on representing the conditions verified by the security chain as a set of logical constraints containing variables. This model is received as input by an SMT solver that checks its satisfiability by searching an assignment of the variables that makes all the constraints true. The second category is based on model checking. In that case, the rules are represented as a finite state machine, and the properties are specified in a (temporal) logic, similar to the verification of the control plane. Starting from this representation, a model checker builds the set of reachable states of the automaton, and checks the validity of properties, building a counter example when the property is violated.

Figure 2 illustrates the interactions among components in our architecture. The security chains that are specified in Pyretic are sent by the security manager to the Synaptic checker. They may first be checked using the Kinetic extension to identify inconsistencies in the control plane. The Synaptic checker then translates the chain specifications into a formal model, which is interpreted by an SMT solver or a model checker. The validation results are then sent back to the security manager, in order to determine whether the security chain can be deployed in the network. If this is the case, the security manager interacts with the SDN controller, which instantiates the chaining with the support of virtual switches.

C. Verification based on SMT solving

We first explain the verification of security chains based on SMT solving. In that context, the specification of a security chain is translated into the SMTlib input language of SMT solvers. Let x be a packet field, x after be the same field after a modification, y be a value, and p 1 and p 2 be given policies. The translation function of Synaptic, in the following denoted by f , is defined in a quite straightforward way based on the grammar of Pyretic, in three main steps:

1) Translation of the Pyretic elementary rules into atomic propositions, as follows:

• f (identity) → true and f (drop) → false,

• f (match(x = y)) → (= x y),

• f (modify(x = y)) → (= x after y). 2) Translation of composition operators (sequence, parallel, and negation) into Boolean expressions, as follows:

• f (p 1 p 2 ) → f (p 1 ) ∧ f (p 2 ), • f (p 1 + p 2 ) → f (p 1 ) ∨ f (p 2 ), 1 • f (∼ p 1 ) → ¬f (p 1 ).

3) Translation of properties of the security chain into

SMTlib constraints. For a given security chain c and a property p, we want to make sure that the implication c → p is valid. This is equivalent to verifying that c ∧ ¬p is unsatisfiable: the latter expression appears in the generated SMTlib input.

In order to illustrate this process, we consider a simple security chain composed of four security functions, noted F i , with i ∈ {1, 2, 3, 4} corresponding to only firewalls. However, our solution is not limited to a specific security function. The Pyretic specification corresponding to this chain appears in Listing 1. Each security function is described by composing the match rules that are applied to ip addresses and port values. For instance, the function F 3 only accepts packets whose source port is 4000, 5000 or 6000. These different functions are then combined using the sequence, parallel and negation operators to specify the final chain. The chain behavior can be described literally as follows: if the security function F 1 accepts a given packet, it is transmitted to the security function F 2 ; otherwise, it is transferred to the security function F 3 . All packets that are accepted by F 2 or F 3 are finally transmitted to F 4 . Assume that we want to check the following property of the chain: each packet accepted by F 1 , F 2 and F 4 or rejected by F 1 and accepted by F 3 and F 4 must be accepted by the chain. (or

(= srcpt port3) (= srcpt port4) (= srcpt port5)) (or (= dstpt port6) (= dstpt port7) (= dstpt port8)))) (not allowed)))) (check-sat)(get-model)(exit)
Listing 2. SMTlib input generated for the toy example (excerpt).

Listing 2 shows the (slightly abridged) result of translating the toy security chain and the above property into the SMTlib language. It includes the declaration of variables and values (abstracting the concrete ip and port values), the translation of the security chain, and the translation of the (negated) property. In particular, each composition operator is translated into a boolean operator, and the overall constraint represents the condition for a packet to be accepted by the chain. The generated file is fed to an SMT solver, which will return a verdict of unsatisfiability if the property is valid or else a model that encodes a counter-example to the property.

D. Verification based on model checking

Our Synaptic checker also supports the verification of security chains based on model checking. This is a direct complement to the existing functionality of Kinetic for verifying the control plane. However, the translation of a Pyretic specification into a finite state machine (or automaton), that can be verified by model checking, is less intuitive than the translation into an SMT model. In order to perform this translation, we extract strictly sequential sub-chains from the security chain described in Pyretic. The following inductive definition introduces the conditions for a chain C to be strictly sequential:

• C ∈ {identity, drop, match, modify}, • C = ∼ C 1 where C 1 is strictly sequential, • C = C 1 C 2
where C 1 and C 2 are strictly sequential.

Subchains built solely from basic rules, negation, and sequential composition are grouped together and define the transitions of the automaton, whereas the parallel composition forms the basis for defining the states of the automaton. Our translation algorithm therefore generates the data plane automaton from the Pyretic specification, based on the following steps:

1) Initially, generate the initial and final states of the data plane automaton. 2) Generate a variable corresponding to each packet attribute (e.g. srcip, srcport) that appears in the specification. These variables are used to specify transition conditions. 3) For each parallel composition, generate an automaton state. These states describe the position of a packet in the security chain. In particular, a packet is considered as accepted when it reaches the final state of the automaton. 4) For each strictly sequential sub-chain, generate an automaton transition. The conditions on this transition correspond to a combination of the elementary rules of the considered sequence.

Figure 3 illustrates the automaton obtained from the Pyretic specification for the example described in the previous subsection. The values for each packet attribute are abstracted by symbolic values, and the states associated to the security chain are noted S 0 , . . . S 5 .

The resulting automaton is represented in a nuXmv model as a finite state machine. It is complemented by the specification of properties to be verified on the chain, written as temporal logic properties to be verified by the model checker. Listing 3 illustrates the properties obtained for our example.

AG(((srcip=ip0

| srcip=ip1) & (srcpt=port0 | srcpt=port1 | srcpt=port2) & (dstpt=port6 | dstpt=port7 | dstpt=port8)) -> EF state=S5) AG(((!(srcip=ip0 | srcip=ip1)) & (srcpt=port3 | srcpt=port4 | srcpt=port5) & (dstpt=port6 | dstpt=port7 | dstpt=port8)) -> EF state=S5)
Listing 3. Generated formal specification interpretable by a model checker.

These properties are expressed in the temporal logic CTL. The formulas are built up using the path quantifiers A and E and the temporal operators F and G. They assert that any data packet whose attributes satisfy the conditions of one of the subchains may reach the final state and therefore be accepted by the security chain. The state machine, together with the temporal formulas, can be passed to nuXmv in order to verify the same properties as the ones verified by SMT solving in the previous sections. Alternatively, properties could be expressed in LTL (linear-time temporal logic), or they could express safety properties that the automaton should satisfy in every state.

E. Implementation of Synaptic

We have implemented our Synaptic checker and its translation algorithms, as a verification package written in Python. It has been designed as an extension of the Pyretic SDN programming language, and exploits the Kinetic extension, part of the Frenetic language family. Our package can be used to verify both the control and data planes of security chains built using the Pyretic language. As depicted in Figure 2, it takes as input the specification of a given security chain expressed in Pyretic, and provides back to the security manager the verification results, corresponding to a verdict on the validity of the chain. Based on this analysis, the security manager can then interact with the SDN controller, so that the security chain is deployed in the programmable network.

Internally, the Synaptic checker relies on two main building blocks: a semantic analyzer capable of interpreting the security chain and its symbols, and a model generator capable of generating a given formal model. The model generator is defined in an abstract manner, so that our checker can easily be extended to other formal languages. We have currently implemented two model generators corresponding to each of the translation algorithms presented below:

• an SMTlib model generator, which produces a formal specification interpretable by an SMT solver, including CVC4 and veriT, • a nuXmv model generator, which produces a formal specification interpretable by the nuXmv solver. 2Our checker can first verify the control plane associated to the security chain, with the support of the Kinetic extension. It then generates the formal specifications corresponding to the data plane, and transmits them to the corresponding solver. The verification of a given security chain can be performed pro-actively by exploring all the possible states reachable by the control automaton, and then all the possible data planes, or can also be performed more reactively by only considering a subset of possible data planes depending on the changes operated on the security chain.

IV. PERFORMANCE EVALUATION

We evaluated the performance of our prototype through an extensive series of experiments. In particular, we wanted to compare the performances obtained with the two translations to SMTlib and nuXmv, as well as for different SMT solvers in order to evaluate the overhead introduced by Synaptic. The experimental setup was based on a MacBook Air laptop computer with a Intel Core i5 (1.7 GHz) processor, and 4Gb of RAM memory. We considered the three following solvers: CVC4 (version 1.4), veriT (version 201506), and nuXmv (version 1.0.1). In order to perform these experiments, we have implemented an additional Python module capable of generating synthetically the security chains that are used as inputs of our Synaptic checker. This generation takes into account several parameters, including:

• the size of the control plane automaton specifying the changes to the security chain in response to network events, measured as the number of states; • the size of the security chain expressed in terms of both width and length; • the number of properties that have to be verified by the checker on the security chain.

While varying these different parameters, we evaluated the response time and memory consumption observed with SMT solving and model checking. This evaluation includes the translation of the specification into a given formal model, and its checking by a given solver. We used the Python time module and the valgrind memory profiler to perform our measurements and obtain statistically grounded results.

A. Impact of the control automaton size

In a first series of experiments, we are interested in quantifying the impact of the size of the control automaton on the verification performance. Our checker supports both the verification of the control and data planes. The size of the automaton directly influences the number of processes required for verifying the security chain. Consider the case of a security chain with an automaton composed of n states. The Synaptic checker therefore must perform n +1 verification processes: one supported by the Kinetic extension to check the control plane, and n processes corresponding to the different automaton states to check the data plane. We use our Python module to generate security chains with different control automaton sizes, varying from 0 to 100 states, in order to quantify the response time and memory consumption induced by our Synaptic checker.

Figure 4 represents the response time (in seconds) of our checker with different backend solvers: CVC4, veriT, and nuXmv. As the verification of the control plane is also based on model checking, it is possible to execute the verification of both control and data planes in a single nuXmv instance. This case corresponds to the last curve, noted nuXmv * . According to these results, we can observe that the response time grows linearly with the size of the control automaton for each of the verification methods. The best performance with that respect is obtained with the nuXmv * configuration, where the control and data planes are verified in the same nuXmv instance. Surprisingly, the worst case is given by the nuXmv curve, with one verification process per automaton state. In that case, the response time is nevertheless still quite acceptable, with a total time of 3.5 seconds for a control automaton of 100 states.

We also quantify the impact on the memory consumption with these different automaton sizes. Figure 5 corresponding to the performance of CVC4, veriT and nuXmv, which are identical. For these three last approaches, the highest memory consumption is always generated by the process supporting the verification of the control plane, using the Kinetic extension and the nuXmv model checker. In these experiments, the nuXmv * approach consumes more memory, with a supralinear behavior when varying the size of the control automaton. In particular, we obtain a memory consumption of 62 Mb for a 100-state automaton. The three other curves (CVC4, veriT, nuXmv) are characterized by a linear behavior, with a memory consumption of about a third of that of nuXmv * for the same 100-state automaton. Therefore, the low response time given by nuXmv * is balanced by a higher memory footprint, which can constitute a limiting factor in scenarios where many security chains have to be checked: for instance, in the case of several security chains protecting several individual applications on a set of smart devices.

B. Impact of the width and length of the security chain

In a second series of experiments, we evaluate to what extent the width and length of the security chain impact the performance of our Synaptic checker. As we are only focusing on the data plane in these experiments, we do not distinguish the nuXmv and nuXmv * cases. We therefore only consider the verification results obtained with the backends CVC4, veriT, and nuXmv.

We first study the impact of the width of the security chain, corresponding to the number of functions (or rules) to be composed in parallel. Figure 6 illustrates the response time of our prototype, while varying the width of the security chain from 100 to 1000. The observed response times grow slightly more than linearly for the three verification methods. We expected the same phenomenon as for the experiments done with the control automaton. In fact, the best performances are obtained with nuXmv. This latter provides a value of 0.8 seconds for a security chain width of 1000, while CVC4 and veriT generate values of respectively 6.4 and 15.9 seconds for the same conditions. When we analyze the results with respect to memory consumption, which are shown in Figure 7, we can observe that veriT requires 62.81 Gb of memory for the widest security chains. CVC4 and nuXmv are much more efficient and only consume respectively 29.89 Mb and 202.67 Mb in the worst cases.

We then focus on the length of the security chain corresponding to the number of functions (or rules) that are composed sequentially in the chain. Figures 8 and9 illustrate respectively the performance in terms of response time and memory consumption, considering a security chain with a length varying from 10 to 100 functions. This parameter has an important impact on the response time of CVC4, with a value of more than 100 seconds in the worst case, while the two other approaches, nuXmv and veriT, maintain acceptable values in these conditions. The results on memory consumption are only indicated for the backends nuXmv and veriT, with CVC4 being disqualified by its bad timing behavior. It appears that veriT consumes less memory than nuXmv for the smallest security chains, from a size of 10 to 60. However, when the length increases from around 80 to 100, nuXmv shows the best results, with a stable value around 26 Mb, while veriT requires almost 127 Mb in the worst conditions. These series of experiments on length and width go in favor of nuXmv (model checking), rather than CVC4 and veriT (SMT solving), when addressing complex security chains.

To complement these results, we perform a third series of experiments, where we quantify simultaneously the impact of the width and length of the security chains for the nuXmv backend, which provided the best performances so far. Figure 10 gives the observed response times with a width varying from 100 to 1000, and a length varying from 10 to 30 (corresponding to the different plotted curves), representing a total number of rules from 1000 to 30000. The most complex security chains require more than 450 seconds of verification time. Some additional optimizations with respect to nuXmv formal modelling could be envisioned to reduce these values, such as reducing the domains associated to packet field representations, or reducing the number of initial states based on the requirements induced by the properties that have to be checked.

C. Impact of the number of properties to be checked

In a last series of experiments, we study more specifically the impact of the number of properties that have to be satisfied by the security chain. Again, due to the previous performance results we only consider the nuXmv backend, where properties are expressed in the temporal logic CTL.

Figure 11 shows both the response time and memory consumption generated by our Synaptic checker with the nuXmv backend, considering a security chain composed of 1000 functions, and varying the number of properties to be checked from 100 to 1000 properties. The performances appear to be linear with respect to the number of properties. More precisely, the response time varies from 0.69 to 1.46 seconds, while the memory consumption goes from 43.34 to 173.99 Mb.

These performance results could be improved depending on the nature of the properties. For instance, when checking simple invariants rather than more complex temporal properties, we could use more efficient dedicated algorithms available in nuXmv.

V. CONCLUSIONS

We have proposed in this paper an automated approach for verifying security chains that are deployed in software-defined networks. Our solution takes into account both the control and data planes of these critical chains, which combine different security functions, such as firewalls, intrusion detection systems, and services for preventing data leakage. We have described the architecture supporting our technique, as well as the behavior and interactions among its different components.

In particular, we have presented our Synaptic checker, which is capable of generating formal models from the security chain specifications, through an extension of the Frenetic language family. We have designed and implemented translation algorithms that underlie verifications of the data plane based on SMT solving and model checking, complementary to the verification of the control plane performed using Kinetic. We developed a prototype of our solution implementing these algorithms, and evaluated its performances through extensive series of experiments based on the backend solvers CVC4, veriT, and nuXmv. The experiments showed the benefits and limits of these methods in terms of response time and memory consumption while varying different sizes of security chains, and numbers of properties to be checked on them.

As future work, we are working on an extension of our translation algorithms to support more complex and advanced rules associated to the security functions. We are also interested in pursuing the optimizations of formal models that we are generating automatically, based on the nature of checked properties. Finally, we want to investigate further the integration of this automated verification process into the context of an automated management framework for security chains.

Figure 1 .

 1 Figure 1. Our architecture supporting the orchestration and verification of security chains.

Figure 2 .

 2 Figure 2. Interaction diagram among the components of our architecture.

F1Listing 1 .

 1 = match(srcip=IP("198.122.37.15")) + match(srcip=IP("253.182.3.14")) F2 = match(srcport=1000) + match(srcport=2000) + match(srcport=3000) F3 = match(srcport=4000) + match(srcport=5000) + match(srcport=6000) F4 = match(dstport=7000) + match(dstport=8000) Pyretic specification of a toy security chain.

(

  set-option:produce-models true) (set-logic QF_LIA) ; Declaration of variables and values (declare-const allowed Bool) (declare-const srcip Int) (declare-const ip0 Int) ... ; Translation of the security chain (assert (and (distinct port3 port5 port0 port7 ip1 port8 port1 port4 port6 port2 ip0) (= allowed (and (or (and (or (= srcip ip0) (= srcip ip1)) (or (= srcpt port0) (= srcpt port1) (= srcpt port2))) (and (not(or (= srcip ip0) (= srcip ip1))) (or (= srcpt port3) (= srcpt port4) (= srcpt port5)))) (or (= dstpt port6) (= dstpt port7) (= dstpt port8)))) ; Translation of the property (and (or (and (or (= srcip ip0) (= srcip ip1)) (or (= srcpt port0) (= srcpt port1) (= srcpt port2)) (or (= dstpt port6) (= dstpt port7) (= dstpt port8))) (and (not(or (= srcip ip0) (= srcip ip1)))

Figure 3 .

 3 Figure 3. Data plane automaton for the toy example.

Figure 4 .

 4 Figure 4. Response time vs. size of the control automaton.

  indicates the maximal memory consumption required during the verification process. Only two curves are visible in this figure, one corresponding to the nuXmv * verification, and another one

Figure 5 .

 5 Figure 5. Memory consumption vs. size of the control automaton.

Figure 6 .

 6 Figure 6. Response time vs. width of the security chain.

Figure 7 .

 7 Figure 7. Memory consumption vs. width of the security chain.

Figure 8 .

 8 Figure 8. Response time vs. length of the security chain.

Figure 9 .

 9 Figure 9. Memory consumption vs. length of the security chain.

Figure 10 .

 10 Figure 10. Response times of nuXmv vs. both width and length.

Figure 11 .

 11 Figure 11. Performance of nuXmv vs. number of properties to be checked.

The output of p 1 + p

is composed of the union of p 1 and p 2 , which explains the translation of this operator as a disjunction.

We use nuXmv rather than its predecessor NuSMV used by Kinetic in order to directly express the constraints on packet attributes.