Aria Adli

Eric Engel

Laurent Romary

Fahime Same

A stand-off XML-TEI representation of reference annotation

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

A stand-off XML-TEI representation of reference annotation

Aria Adli 1 , Eric Engel 1 , Laurent Romary 2 , Fahime Same 1 1 University of Cologne, 2 INRIA {aria.adli, eric.engel, f.same}@uni-koeln.de, laurent.romary@inria.fr What's new?

.

. § Separation of primary data (transcription) and secondary data (annotations) § Incremental build-up of annotation layers § Three-way distinction of reference-related features § Local grammatical features (grammatical gender, number, person) and sentence-bound discourse features (information status) apply to referring expressions. § Inherent semantic features of referents (animacy, cardinality, specificity, concreteness/abstractness) apply to the referents.

Their values are constant over the entire discourse. § Relations between referents (subset, partOf) are recognized as permanent, allowing for a better representation of so-called bridging phenomena in discourse.

Model features

Background and previous work § User-friendly visualization and query interface (in ANNIS) § Integration of different levels of information structure analysis § [±topical] as feature of a referring expression? § Markup of focus domains directly on syntactic phrases? § Further expansion of the multilayer architecture to the analysis of other cross-sentential phenomena, e.g., rhetorical structure, dialogue analysis

  [1] ISO 24611:2012 Language resource management -Morpho-syntactic annotation framework (MAF). [2] ISO 24615-1:2014 Language resource management -Syntactic annotation framework (SynAF) -Part 1: Syntactic model. [3] ISO/AWI 24617-9 [Under development] Language resource management -Semantic annotation framework -Part 9: Reference Annotation Framework. [4] Poesio, M., 2004. The MATE/GNOME proposals for anaphoric annotation, revisited. In Proceedings of the 5th SIGdial Workshop on Discourse and Dialogue at HLT-NAACL 2004, 154-162. [5] Riester, A. and Baumann, S., 2017. The RefLex scheme -Annotation guidelines. SinSpeC 14. 1-31. [6] Salmon-Alt, S. and Romary, L. 2005. The Reference Annotation Framework: A case for semantic content representation. In H. Bunt (ed.), IWCS-6, Tilburg, Netherlands: ACL SIGSEM. References § Reference annotation is an essential element in corpus-linguistic research at the interface of syntax and discourse structure. § We propose a representation based on a semantically plausible model distinguishing referring expressions, discourse referents, and links. § Neater representation of grammatical vs. semantic features and of phenomena of indirect reference ("bridging") § More efficient annotation of longer stretches of natural discourse § Integrated in a multilayer stand-off annotation in compliance with TEI-ISO standards: MAF [1], SynAF [2], RAF [3]

  Multilingual corpus project on French, Spanish, Persian, and Catalan § Spontaneous speech data elicited in a game task § More info: www.sgscorpus.com § Building on extensive work on (stand-off) representation of referring expressions and referential links [4,6] and types of links [5] § So far, no existing models making a clean distinction between expression-related and entity-related features and links Example


L Gremillard
email: laurent.gremillard@insa-lyon.fr

K Biotteau-Deheuvels

P Clement

M Perron

P Duvauchelle

S Meille

J Chevalier

Aatt

Keywords: thermal treatment, zirconia-toughened alumina, simulation, sinter/sintering

During the sintering process of ceramic pieces, temperature gradients may arise inside the piece if the thermal cycle is too fast with regards to the size of the ceramic body. Too large temperature gradients may lead to density gradients, thus to strain gradients which may result in defects or in the fracture of the pieces. This article presents both an experimental setup and the first steps towards a numerical method to evaluate these thermal gradients and their evolutions during sintering, in the case of large zirconia, alumina and zirconia-toughened alumina balls. Finally it discusses the advantages of Rate Controlled Sintering thermal cycles in terms of diminishing the temperature gradients during densification.

Introduction

The aim of Total Hip Arthroplasty is to rebuild a functional hip joint. Charnley first used a femoral head with a diameter of 41.5 mm [START_REF] Wroblewski | The Charnley Hip Replacement -43 Years of Clinical Success[END_REF]. At that time (the early 1960's) he unfortunately had to abandon this large head size because of the limitations of the materials wear performances. After a trend to decrease the size of the femoral heads (down to 22.2 mm diameter, allowing lower wear debris generation) larger and larger femoral heads have been recently developed. They now reach 46 mm diameter (sometimes up to 51), which allows more flexibility for the surgeon, more comfort for the patient and lower risks of dislocation. New materials and techniques serve this will of increasing the femoral head diameter toward the size of the natural hip joint (usually between 45 and 55 mm). The use of ceramic-onceramic bearing couples give most promising results compared to other bearing couples (lower friction, very few wear debris). Moreover the third generation ceramic-on-ceramic articulation could be of particular benefit for young and/or active adults, with the use of large ceramic femoral heads [START_REF] Byun | Third-Generation Ceramic-On-Ceramic Total Hip Arthroplasty in Patients Younger Than 30 Years with Osteonecrosis of Femoral Head[END_REF]. This diameter evolution is supported by many studies insisting on the improvement of the living quality: increased Range Of Motion [START_REF] Burroughs | Range of Motion and Stability in Total Hip Arthroplasty With 28-, 32-, 38-, and 44-mm Femoral Head Sizes: An In Vitro Study[END_REF], decreased risk of dislocation [START_REF] Burroughs | Range of Motion and Stability in Total Hip Arthroplasty With 28-, 32-, 38-, and 44-mm Femoral Head Sizes: An In Vitro Study[END_REF][START_REF] Jameson | Lower rates of dislocation with increased femoral head size after primary total hip replacement: a five-year analysis of NHS patients in England[END_REF][START_REF] Triclot | Update -"Big-head": The solution to the problem of hip implant dislocation?[END_REF][START_REF] Conroy | Risk Factors for Revision for Early Dislocation in Total Hip Arthroplasty[END_REF], decrease of wear debris amount [START_REF] Slonaker | Review of wear mechanisms in hip implants: Paper IIceramics IG004712[END_REF], increased stability of the hip [START_REF] Burroughs | Range of Motion and Stability in Total Hip Arthroplasty With 28-, 32-, 38-, and 44-mm Femoral Head Sizes: An In Vitro Study[END_REF] and decreased risk of impingement [START_REF] Kluess | Influence of femoral head size on impingement, dislocation and stress distribution in total hip replacement[END_REF][START_REF] Marchetti | Component impingement in total hip arthroplasty: Frequency and risk factors. A continuous retrieval analysis series of 416Â cup[END_REF]. However large technical ceramic pieces face an important, intrinsic processing difficulty: their large size and low thermal conductivity make them susceptible to the apparition of potentially large thermal gradients during sintering. Moreover, industrial ceramic components are commonly produced by Cold Isostatic Pressing (CIP) or Uniaxial Pressing, followed by sintering in air. This process can induce defects such as density gradients in the component

Gremillard et al. Thermal gradients during sintering of ceramic balls 2

which will be exaggerated by sintering [START_REF] Green | Constrained sintering: A delicate balance of scales[END_REF]. In fact, sintering leads to consolidation of the material through particle bonding and porosity elimination. During the heat treatment, mass transfer takes place to minimize surface energies and large microstructural changes occur, leading to densification and grain growth. Since temperature directly influences the shrinkage rate, the generated thermal gradients linked to the relatively low density gradient may result in differential shrinkage throughout the piece, thus distortion, potential cracking or even premature end of the shrinkage [START_REF] Missiaen | Modelling of sintering: recent developments and perspectives[END_REF]. This is even more pronounced with the application of recent, fast sintering techniques such as Spark Plasma Sintering and Microwave Sintering. Consequently, successful near net shape forming of ceramic powder based components can be made easier by a more careful consideration of thermal gradients, through both experimental and simulation works.

Since the 1960's, many researches have been conducted to describe and/or model the sintering mechanisms and processes. Modeling of the sintering process is thoroughly investigated, via an important numbers of methods at different scales (mostly scale of the grains, and scale of the parts). The most advanced studies take into account multi-scale modeling [START_REF] Green | Constrained sintering: A delicate balance of scales[END_REF][START_REF] Missiaen | Modelling of sintering: recent developments and perspectives[END_REF][START_REF] Maximenko | Direct Multi-Scale Modeling of Sintering[END_REF]. Macroscopic models of sintering generally aim to describe shape evolution and mechanical constraints at large scale, with the hypothesis of a homogeneous microstructure by finite elements methods (FEM). More recent developments also take into account the heterogeneity of the initial microstructure. For example, X-Ray Computed Micro-Tomography (XRCT) [START_REF] Bernard | First direct 3D visualisation of microstructural evolutions during sintering through X-ray computed microtomography[END_REF] combined with modeling revealed itself a useful tool to obtain 3D visualization of microstructural evolutions during isothermal sintering. These observations served as a basis for the analysis of sintering phenomena in order to extract more statistical results and experimental data to improve modeling. [START_REF] Chaix | Quantitative analysis of microstructure and modeling of sintering[END_REF] As an example, Schoenberg et al. [15,[START_REF] Schoenberg | Stresses and distorsion due to green density gradients during densification[END_REF] recently used finite element analysis to calculate the stresses resulting from the variations in green density. The approach involved measurements of the viscosity constants as a function of green density, combined with X-ray computed tomography to describe the density variations prior to sintering. The finite element analysis then allowed a simulation of the stresses and strains occurring in the densification process of a specific sample. This allowed regions of high stress to be identified and this was correlated to the observed damage.

However, only few of these studies even mention the existence of a thermal gradient during sintering [START_REF] Vallauri | Simulation of Solid State Sintering through FE Modeling for the Optimum Design of 3D Parts[END_REF][START_REF] Zhang | Stress evolution in graded materials during densification by sintering processes[END_REF]. The goal of this paper is thus to characterize experimentally the thermal gradients in large ceramic spheres during sintering and make a first attempt at simulating them using a finite element model. The latest would give information about the global shrinkage kinetics. This is a first step towards a finite element simulation that could also predict distortions of the compact during and after sintering and internal stresses during sintering.

Materials and methods

Finite Element Model

Finite Element Simulations (with the software ANSYS 14.5) of the heating and plateau phases of a sintering thermal cycle were conducted on 51 mm diameter spheres (diameter before simulation) made of zirconia-toughened alumina, monolithic alumina or zirconia (3Y-TZP) ceramics. The calculations were conducted on an axisymmetric 2D model with a regular mesh (PLANE13 quadrangle elements [START_REF]ANSYS, 4.13 PLANE13 2-D Coupled-Field Solid (UP19980821[END_REF], (2-D coupled field solid)), coupling mechanical and thermal analyses in transitory regime. The output of the calculation was the temperature field in the ceramic sphere versus time. From this output, the evolutions with time of three parameters were specially considered: the temperature on the center of the sphere (T C ), on its surface (T S ) and the difference between those (T=T S -T C ). Since we were only interested in Gremillard et al.

Thermal gradients during sintering of ceramic balls 3

the temperature gradient inside the sphere, the issues linked to heat transfer between the furnace and the sphere were ignored using Dirichlet conditions, i.e. by setting the sphere surface temperature to the furnace temperature. All materials properties necessary for the simulation vary during sintering, both with temperature and density. They are represented on Figure 1, and explained below. The properties of ZTA were approximated using a simple proportionality between those of zirconia (10%) and those of alumina (90%). The evolution of density versus temperature, (T), was determined experimentally by dilatometry on small cylinders ( 1 cm diameter, 5 mm thickness), following a typical thermal sintering cycle (1K•min -1 until 473 K, 2h dwell, 5 K•min -1 up to 1793 K, 2 h dwell and ~5 K•min -1 down to room temperature), giving directly access to the variations of porosity with temperature P(T). It was then translated in the simulation by assimilating the density variation with temperatures to a thermal shrinkage. A secant apparent thermal expansion coefficient varying with temperature, () T  , was therefore defined according to eq. 1:

          00 00 L T L T T T T 3 T T        eq.1
Where L(T) is the length of the sample used for the dilatometry measurement at a temperature T, L 0 is its initial length (at room temperature T 0 ), and  and  0 are the density of the same sample at temperature T and room temperature. Eq.1 is valid for an isotropic shrinkage, which is reasonable for samples pressed by CIP. The increase in density with temperature results in a negative apparent thermal expansion coefficient. One should note that the densification versus temperature curve is strongly dependent on the powder and pressing parameters, and may vary greatly from one powder to another even if they have the same nominal composition. Conversely, powders of different nature may present the same densification versus temperature curve. Therefore, for the sake of comparison, the same densification vs temperature curve was used for all materials considered here.

Gremillard et al. Thermal gradients during sintering of ceramic balls 4

All other material parameters were expressed as a function of density (or porosity) and temperature (when applicable). As density itself is expressed as a function of temperature, the variation of every single parameter during simulation can be expressed as a function of temperature only (but this variation accounts for both the effects of temperature and density). This approach greatly simplifies the simulation inputs.

The material was considered as purely linear-elastic. Although this is a huge approximation, especially at higher temperatures, this should not influence the thermal gradients. Young's modulus was considered as varying linearly with porosity, following Dewey-Mackenzie relations [START_REF] Pabst | Elasticity of porous ceramics-A critical study of modulus-porosity relations[END_REF] as:

         00 00 E 3 1 9 T T E 5 1P 2 7 5    νν ν eq. 2
Where E(T) is the Young's modulus during the sintering cycle, E 0 the Young's modulus of the 100% dense material at T 0 ,  0 =0.26 the Poisson's coefficient, and P(T) the porosity of the material (varying with temperature). E 0 was considered as constant with temperature [START_REF] Pabst | Elasticity of porous ceramics-A critical study of modulus-porosity relations[END_REF]: 400 GPa for alumina et 200 GPa for zirconia. Although in truth E 0 varies slightly with temperature [START_REF] Mignard | Thermal shock behaviour of a coarse grain porous alumina[END_REF], this variation does not impact the thermal behavior.

The thermal conductivity

K(T) (W•m -1 •K -1
) is described by equation 3(a), [START_REF] Litovsky | Gas Pressure and Temperature Dependences of Thermal Conductivity of Porous Ceramic Materials: Part 1, Refractories and Ceramics with Porosity below 30%[END_REF] in which K D (T) is the thermal conductivity of the dense material (depends on temperature), and n is a dimensionless coefficient chosen equal to 8 [START_REF] Rhee | Porosity -Thermal conductivity correlations for ceramic materials[END_REF]. K D varies with temperature following the phenomenological equations 3(b) [START_REF] Auerkari | Mechanical and Physical Properties of Engineering Alumina Ceramics[END_REF], and 3(c) [START_REF] Auerkari | Mechanical and Physical Properties of Engineering Alumina Ceramics[END_REF][START_REF] Swain | Thermal diffusivity, heat capacity and thermal conductivity of porous partially stabilized zirconia[END_REF][START_REF] Mévrel | Thermal diffusivity and conductivity of Zr1-xYxO2-x/2[END_REF][START_REF] Limarga | The grain size and temperature dependence of the thermal conductivity of polycrystalline, tetragonal yttria-stabilized zirconia[END_REF] for alumina and zirconia respectively (T expressed in Kelvin). These equations are derived from the data measured by the cited authors. Note that K(T) varies both with porosity and temperature. However, since a dependence of porosity on temperature only was postulated, K(T) can be expressed (in this simulation) as varying only with temperature.

        D 2 1P KT 1 n P T K T T     eq. 3(a)       3 D K T 5.5 34.5exp 3.3 10 T 273       eq. 3(b)       3 D K T 1.89 0.974exp 1.85 10 T 273       eq. 3(c)
The specific heat capacities C P of alumina and zirconia were respectively calculated using the approximations shown in eq. 4(a) [START_REF] Miyayama | Engineering properties of single oxides[END_REF] and 4(b) [START_REF] Swain | Thermal diffusivity, heat capacity and thermal conductivity of porous partially stabilized zirconia[END_REF] (T in Kelvin and C p in J•kg -

     eq. 4(b)
The properties of zirconia-toughened alumina composites were calculated using a simple rule of mixtures.

Different thermal cycles were used during the simulations, in order to evaluate their effect on temperature homogeneity in the spheres. The first one (EC, standing for Experimental Cycle) reproduces the sintering cycle used for the experiments. Others reproduce typical sintering cycles, with heating rates of [START_REF] Triclot | Update -"Big-head": The solution to the problem of hip implant dislocation?[END_REF][START_REF] Green | Constrained sintering: A delicate balance of scales[END_REF][START_REF] Pabst | Elasticity of porous ceramics-A critical study of modulus-porosity relations[END_REF]30,40 or 50 K•min -1 until 1793 K. The last two reproduce Rate Controlled Sintering cycles (5RCS and 19RCS: respectively 5 or 19 K•min -1 up to 1430 K, followed by a heating rate proportional to the densification rate so that the duration of the 19 RCS cycle is the same as that of a typical 5 K•min -1 cycle).

A convergence analysis was conducted on a sintering cycle at 10 K•min -1 applied to zirconia, with a number of elements between 700 and 7000. In these calculations the maximum difference between the temperatures on the surface and in the center of the sphere (T Max ) was picked up. Figure 2 shows a maximum difference of 0.5 K between the different T Max calculated with the different meshes. Thus it was decided to use a 1100 elements mesh, providing both adequate accuracy and fast simulations. slow ramp (1K•min -1 ) followed by a dwell at 473 K, then a 5 K•min -1 up to 1793 K, followed by a 2 h dwell and a ~5 K•min -1 cooling ramp, which is consistent with sintering cycles used in real industrial conditions. The experimental setup is described on Figure 3, together with the applied thermal cycle.

Results

Finite Element Model

The temperature differences between the surface and the center of a sphere (T=Ts-Tc) obtained with the FE simulations are presented in Figure 4 for alumina and zirconia (simulations for ZTA do not differ significantly from those of alumina). As expected from their relative thermal conductivity, for all thermal cycles the temperature inside zirconia spheres is more heterogeneous than inside alumina or ZTA spheres. The simulations with Rate-Controlled Sintering conditions show that at the beginning, during the fast heating rate T logically increases at the same pace as it would with a 5K•min -1 or 19K•min -1 heating rate. However, when densification starts and thus the heating rate decreases (down to near 1 K•min -1 ), T dramatically decreases and becomes the lowest of all simulated conditions: in average during the densification steps (between 1440 K and 1793 K),

T is 68% lower during Rate Controlled Sintering (whatever the initial heating rate) than during a sintering at 5 K•min -1 constant heating rate for alumina, and 64% lower for zirconia (Figure 6). A comparison of T obtained for a 5K•min -1 cycle and a 19RCS cycle in zirconia between 100 and 1700 K is shown on figures 7, Reducing the initial heating rate (5RCS cycle) further reduces the temperature gradient before densification starts, without affecting it after. This is illustrated on zirconia on figure 8.

Gremillard et al.

Thermal gradients during sintering of ceramic balls 8

In-situ measurement of temperatures during sintering

The temperature differences between the surface (T S ) and center (T c ) of the spheres were only measured on ZTA. Erreur ! Source du renvoi introuvable.The results are presented on Figure 9. A small difference exists at low temperature, but is removed by the plateau at 473 K. After this plateau, the temperature difference increases steadily until it reaches almost 50 K at an external temperature around 800 K. Then it mostly decreases until it reaches 30 K at an external temperature of 1600 K. Finally, the temperature difference decreases in two significant steps: the first step takes place at an external temperature of 1650 K (internal temperature of 1620 K), when a significant densification has already taken place in the whole body (porosity below 30 %) and thus the thermal conductivity increases (arrow indicated in Figure 9). The second decrease takes place when the external temperature reaches the soaking temperature (1793 K) and stops increasing, and the temperature naturally homogenizes with time in the body.

Gremillard et al.

Thermal gradients during sintering of ceramic balls 10

Discussion

In general, the orders of magnitude of the temperature differences between the surface and center of the spheres computed using the finite element model are in good agreement with the measured ones, for similar conditions. Given the difficulty of the experiments, it is difficult to know whether the experiments or the model are most precise. Indeed, the experimental measurement may be impaired by many factors. These factors are:

-The electromagnetic fields inside the furnace that can change the potentials in the thermocouples.

-The thermal conductivity of the thermocouples can also change the measured temperature, mostly in the center of the heads (thus the temperature of the center would be overestimated during heating, and the measured thermal gradient artificially decreased).

-The thermocouples are bent during the CIP step; this can affect their metallurgical state and thus the temperature measurement.

The most important finding of this study is perhaps the demonstration of the existence of very low thermal gradients during densification when using Rate Controlled Sintering thermal cycles. RCS thermal cycles should be separated in two different phases: heating and densification. During the heating phase, the first heating rate can be large (it was limited here to 19 K / min), and give rise to high temperature differences between the surface and the core (44 K in alumina, 115 K in zirconia, both at a surface temperature of 1433 K, for a 19RCS cycle). Since the temperatures during the heating phase are compatible with an elastic behavior of the materials, with the Young's modulus given in Figure 1 (or Eq. 2) (resp. 84 and 42 GPa for alumina and zirconia), and using a Coefficient of Thermal Expansion of 6.3•10 -6 for alumina and 11 •10 -6 for zirconia, one can calculate maximum elastic stresses of about 23 MPa for alumina and 53 MPa for zirconia. Of course these values should be taken only as orders of magnitudes. But they still exceed by far the strength of green ceramic bodies (a few MPa at most [START_REF] Pizette | Green strength of binder-free ceramics[END_REF]). Thus heating at such a high rate may result in cracking of the body. However, the initial heating rate applied here for 19RCS was calculated so that the total duration of RCS thermal cycle is the same as the duration of a 5K/min thermal cycle. Since the initial heating rate has no influence on densification, it is easy to decrease it to avoid this problem, as it is shown by the 5RCS cycles. As soon as the densification phase starts and microstructures start to develop, the thermal gradients reach very low values (maximum 2 K for alumina and 5K for zirconia above 1573 K, when densification takes place) since simultaneously the heating rate becomes low and the conductivity starts to increase. This presents two advantages: it should produce pieces with very low distortions and with very homogeneous microstructures (because of the very low thermal gradients, every location of the body will experiment the same thermal history). Thus an ideal sintering cycle for large ceramic bodies, in terms of homogeneity and distortions, may be constituted by a slow heating rate up to the temperature where densification starts (1560K here) followed by a ratecontrolled sintering cycle during densification. Some other results of the finite element analysis are expected, such as the proportionality of the maximum thermal gradient with the heating rate. However, the (applied) temperature for which the maximal thermal gradient is observed (T(T Max )) lays between 1500 and 1700K, depending on the material and the sintering conditions. This indicates first that the temperature gradient decreases between the temperature at which the densification starts (~1500 K) and the soaking temperature (1793 K). This is explained by an increase of the thermal conductivity due to densification. Densification starts at the surface, where the temperature is the highest, and the densification front "propagates" to the interior. When density is high enough everywhere in the sintering body, the thermal gradient can start to decrease. This is reached for an internal temperature of around 1500 K for the three materials, corresponding to a 2% density increase (from 60 to 61.2% theoretical density in the center of the pieces).

As for now, this purely phenomenological model shows three main limitations: -The simulations do not take into account the elasto-visco-plastic behavior of the materials, which prevents an analysis of internal stresses. However, the model could be improved with further development based on experimental data.

-An important point for a more accurate modeling of the thermal gradient is accounting for density inhomogeneity. Indeed, all the simulations presented here consider a perfectly homogeneous starting material. However, in green bodies, especially in large ones, particles may be heterogeneously packed. These heterogeneities can affect densification either directly or indirectly. The direct effect results from density itself: the less dense parts will sinter with more difficulty, and be submitted to a larger shrinkage (thus larger distortions). The indirect effect is a thermal effect: less dense parts have a lower thermal conductivity and will be submitted to larger thermal gradients, thus to larger distortions.

-The influence of heating rate on the sintering kinetics was neglected. It is known that faster heating rate can result in faster densification (not taking into account the thermal gradients). This effect is neglected here and could be accounted through the establishment of Master Sintering Curves.

Other, less important, points could be more easily refined:

-Most physical characteristics of materials were not measured, but extrapolated from literature.

-The evolution of density with temperature was measured on ZTA only, and was used for the simulations for the three materials. This is of limited importance, since the densification behavior varies greatly with the pre-sintering processing steps (pressing, casting…) and the nature of the powder (particle size, dopants…).

-The simulations are conducted only on the heating part of the cycle and on the plateau, preventing the analysis of thermal gradients during cooling. Simulating the cooling would necessitate to inject the results obtained by the present model after the plateau into a second

Gremillard et al. Thermal gradients during sintering of ceramic balls 12

model in which the density remains constant with temperature (unlike during heating) while all other physical properties vary with temperature. This is also of limited importance, since cooling in industrial furnaces is generally very slow due to the size and thermal inertia of the furnaces.

The main interest of running such simulation is to obtain a useful and effective tool in order to simulate the sintering behavior of large ceramic parts. This model is quite simple and easily transferable to other geometries but also to other materials. Materials properties and sintering conditions can be easily modified to adapt in every situation. Finally the calculation time is relatively fast, as it takes less than one hour on a commercial laptop, which put this model as a useful tool for calculating sintering thermal cycles.

Conclusion

This article aimed at establishing a macroscopic finite element model of sintering, and comparing these simulations to the measured thermal gradients inside a ceramic sphere during sintering. Both experiments and model give results with comparable orders of magnitudes.

The model shows that the temperature gradients inside a large ceramic piece vary linearly with the heating rate, and that the maximum gradient is reached just above the temperature where shrinkage begins. Moreover, it may allow the determination of optimal sintering cycles taking into account the temperature (thus microstructural) homogeneity inside the pieces. It might however not be sufficient in its current version to predict the shape changes during sintering, as these changes are related to temperature gradients via shrinkage and mechanical stresses, which are highly influenced by the thermo-mechanical behavior of the materials (approximated here as purely elastic). ,nbligne,1,1, MPTEMP,,298.15,323.15,358.15,373.15,398.15,423.15 ! Variation de la température MPTEMP,,448.15,473.15,498.15,523.15,548.15,573.15 ! Pour propriétés mécaniques et thermiques MPTEMP,, 598.15, 623.15, 648.15, 673.15, 698.15, 723.15 MPTEMP,, 748.15, 773.15, 798.15,823.15,848.15, 873.15 MPTEMP,, 898.15, 923.15, 958.15, 973.15, 998.15, 1023.15 MPTEMP,, 1048.15, 1073.15, 1098.15,1123.15,1148.15, 1173.15 MPTEMP,, 1198.15, 1223.15, 1248.15, 1273.15, 1298.15, 1323.15 MPTEMP,, 1348.15, 1373.15, 1398.15,1423.15, 1428.15, 1433.15 MPTEMP,, 1438.15,1443.15,1448.15,1453.15,1458.15, 1463.15 MPTEMP,, 1468.15,1473.15,1478.15,1483.15,1488.15, 1493.15 MPTEMP,, 1498.15,1503.15,1508.15,1513.15,1518.15, 1523.15 MPTEMP,, 1528.15,1533.15,1538.15,1543.15,1548.15, 1553.15 MPTEMP,, 1558.15,1563.15,1568.15,1573.15,1598.15, 1623.15 MPTEMP,, 1648.15, 1673.15, 1698.15,1723.15,1748.15, 1773.15 MPTEMP,, 1793 PLVAR,[START_REF] Burroughs | Range of Motion and Stability in Total Hip Arthroplasty With 28-, 32-, 38-, and 44-mm Femoral Head Sizes: An In Vitro Study[END_REF][START_REF] Jameson | Lower rates of dislocation with increased femoral head size after primary total hip replacement: a five-year analysis of NHS patients in England[END_REF]FILLDATA,192,,,,0,0 FILLDATA,193,,,,[START_REF] Wroblewski | The Charnley Hip Replacement -43 Years of Clinical Success[END_REF]0 FILLDATA,194,,,,0 FILLDATA,195,,,,[START_REF] Wroblewski | The Charnley Hip Replacement -43 Years of Clinical Success[END_REF]1 VARNAME,195,NSET

Figure 1 :

 1 Figure 1: Materials properties used for the simulation vs temperature (blue squares: zirconia; black circles: alumina. (a): Specific heat capacity (Cp, J•kg-1•K-1); (b): Thermal conductivity (K, W•m-1•K-1); (c): porosity (measured by dilatometry, the same for all materials); (d): Young's modulus (E, GPa).

Figure 2 :

 2 Figure 2: (a) Convergence analysis giving the maximum temperature difference between the center and the surface for different number of elements (left axis) and the difference in this value as compared to the case of 7000 elements (right axis), and (b) representation of the meshing with 1100 elements.

Figure 3 :

 3 Figure 3: Experimental setup (a) and temperature cycle (b) used to measure the thermal gradients in a sphere during sintering.

Figure 4 :

 4 Figure 4: T as a function of temperature simulated in alumina (left) and zirconia (right).

Figure 5 :

 5 Figure 5: (a) Evolutions of TMax and T(TMax) (temperature for which T is maximum) with heating rate, and (b) proportionality between TMax and T(TMax) (the size of the bubbles is proportional to the heating rate).

Figure 6 :

 6 Figure 6: Comparison between Rate-Controlled Sintering (19RCS, thick lines without symbols) and sintering at 5K•min-1 (thin lines with symbols) for alumina (black) and zirconia (blue). Red lines indicate the temperature

Figure 7 : 9 Figure 8 :

 798 Figure 7: Comparison between Rate-Controlled Sintering at two initial heating rates: 19RCS (thick lines) and 5 RCS (thin lines) for zirconia (blue). Red lines indicate the temperature (thick, dashed line: 19RCS; thin, dashed line: 5RCS)

Figure 9 :

 9 Figure 9: Temperature difference between the surface (Ts) and the center (Tc) of a ZTA sphere measured in situ during densification (open diamonds) and calculated using the finite element model (full diamonds).

15 !

 15 , , *TREAD, tempreelsurfAnsys,'tempreelsurfAnsys','txt',, , !****************************************************************************** ! Variables géométriques et données du problème !******************************************************************************************************************************************************** ! Création de la géométrie et choix du maillage !****************************************************************************La température initiale est de 25°C Gremillard et al. Thermal gradients during sintering of ceramic balls 16 ! ****************************************************************************** ! Propriétés matériau ! ******************************************************************************

  .15 mp,dens,1,D ! Application des propriétés du matériau -densité mp,PRXY,1,PR ! Application des propriétés du matériau -coefficient de Poisson ! Application des propriétés du matériau -module d'Young MPDE,EX,1 ! Publication Latella et Liu AMESH,ALL ! ****************************************************************************** ! Résolution ! ****************************************************************************régime transitoire (dépend du temps) OUTRES,ALL,ALL ! Accessibilité des résultats à chaque incrément de temps KBC,0 ! Application progressive du chargement thermique (rampe) ! Conditions aux limites et symétries nsel,all ! On bloque la translation selon x sur l'axe de symétrie LSEL,la ligne à x=0 nsel,r,loc,y,0 ! Sélection d'un point à y=0 dans le ligne D,ALL, , , , , ,UY ! Le centre du cercle est bloqué aussi selon y nsel,all ! Application de la température sur le contour -première phase LSEL,S, , ,1 NSLL,S,1 D,ALL,TEMP,%tempreelsurfAnsys% ! Application de la température grâce au fichier extérieur nsel,all Gremillard et al. Thermal gradients during sintering of ceramic balls 23

  •K -1 ). It was considered as constant vs. density.

	  C T 1.126 10 +1.256 10 T 31 p    		2 3.476 10 T 	7	eq. 4(a)
	  C T 5.608 10 +1.104 10 T 21 p 		2 1.050 10 T	7

Example of Ansys code used for the simulations

	! Modèle de frittage d'une tête de prothèse de hanche en zircone
	! Etude thermique
	! Cp tient compte de l'évolution avec la température
	! K lié à température et densité ; K liste au lieu de polynôme
	! rampe de montée à 5°C/min au lieu de 19
	nbligne =883	! Lecture du fichier de température
	*DIM, tempreelsurfAnsys,TABLE

Gremillard et al. Thermal gradients during sintering of ceramic balls 15

Annex 1:

Acknowledgements

The authors wish to acknowledge the French ANR (Agence Nationale pour la Recherche) for its support to project OptHip under grant number MAPR07-0014. The authors also thank Medical Group (Vaulx en Velin, France) for providing ZTA powder, and François Valdivieso (Ecole des Mines de St Etienne) for useful discussions around this topic.

MPDE,KXX,1 MPDATA,KXX,1,,0.744205666362131 MPDATA,KXX,1,,0.739095596637678 MPDATA,KXX,1,,0.725317838086449 MPDATA,KXX,1,,0.712389553871171 MPDATA,KXX,1,,0.702729617718701 MPDATA,KXX,1,,0.693507008920525 MPDATA,KXX,1,,0.684701928669299 MPDATA,KXX,1,,0.679542764524676 MPDATA,KXX,1,,0.677697314511143 MPDATA,KXX,1,,0.675818250744359 MPDATA,KXX,1,,0.670356793401028 MPDATA,KXX,1,,0.664776389295863 MPDATA,KXX,1,,0.657917603811226 MPDATA,KXX,1,,0.649814852874495 Gremillard et al. Thermal gradients during sintering of ceramic balls 18 MPDATA,KXX,1,,0.645117513829061 MPDATA,KXX,1,,0.638952091194269 MPDATA,KXX,1,,0.633450167613004 MPDATA,KXX,1,,0.628009591672646 MPDATA,KXX,1,,0.619845852811009 MPDATA,KXX,1,,0.61491038470031 MPDATA,KXX,1,,0.610785359834162 MPDATA,KXX,1,,0.605699652557927 MPDATA,KXX,1,,0.601404615215334 MPDATA,KXX,1,,0.597304024014825 MPDATA,KXX,1,,0.593389075935807 MPDATA,KXX,1,,0.58871248940573 MPDATA,KXX,1,,0.581457272872958 MPDATA,KXX,1,,0.57989661684504 MPDATA,KXX,1,,0.576659431032107 MPDATA,KXX,1,,0.573568800008418 MPDATA,KXX,1,,0.570172038228086 MPDATA,KXX,1,,0.567800963225432 MPDATA,KXX,1,,0.565111375255013 MPDATA,KXX,1,,0.562103812238947 MPDATA,KXX,1,,0.560091978145506 MPDATA,KXX,1,,0.557751393522115 MPDATA,KXX,1,,0.555516772506007 MPDATA,KXX,1,,0.553915890292847 MPDATA,KXX,1,,0.55134644961921 MPDATA,KXX,1,,0.551511950850398 MPDATA,KXX,1,,0.54982190005783 MPDATA,KXX,1,,0.548388382263173 MPDATA,KXX,1,,0.546170019417321 MPDATA,KXX,1,,0.545933669618772 MPDATA,KXX,1,,0.543171231221213 MPDATA,KXX,1,,0.542034790840639 MPDATA,KXX,1,,0.541746973740741 MPDATA,KXX,1,,0.54232522225666 MPDATA,KXX,1,,0.542906440599495 MPDATA,KXX,1,,0.541680029080365 MPDATA,KXX,1,,0.541402283868441 MPDATA,KXX,1,,0.542935839470111 MPDATA,KXX,1,,0.543527201200008 MPDATA,KXX,1,,0.544121463766273 MPDATA,KXX,1,,0.544718615377959 MPDATA,KXX,1,,0.544322254731248 MPDATA,KXX,1,,0.544057451718477 MPDATA,KXX,1,,0.545079619849563 MPDATA,KXX,1,,0.545976347323049 MPDATA,KXX,1,,0.548328383999059 MPDATA,KXX,1,,0.550691586307503 MPDATA,KXX,1,,0.55120672856697 MPDATA,KXX,1,,0.552810715437804 MPDATA,KXX,1,,0.557553408431787 MPDATA,KXX,1,,0.561127200708515 MPDATA,KXX,1,,0.564742636711708 MPDATA,KXX,1,,0.568550214121921 MPDATA,KXX,1,,0.572363202535151 MPDATA,KXX,1,,0.582765837544311 MPDATA,KXX,1,,0.590612243219992 MPDATA,KXX,1,,0.595737540998404 MPDATA,KXX,1,,0.600980012767243 MPDATA,KXX,1,,0.618090106629032 MPDATA,KXX,1,,0.625828480050249 MPDATA,KXX,1,,0.642760368517748 MPDATA,KXX,1,,0.674275541997735 MPDATA,KXX,1,,0.756923158724945 MPDATA,KXX,1,,0.89831633351243


Abdelouahab Attia
email: attia.abdelouahab@gmail.com

Abdelouahab Moussaoui

Youssef Chahir

An EEG-fMRI Fusion Analysis Based on Symmetric

Keywords: EEG-fMRI data fusion, symmetric as-symmetric approach, BOLD, HRF, Dempster Shafer theory, PHA method

  L'archive ouverte pluridisciplinaire

Introduction

Electroencephalogram (EEG) is a non-invasive technique that measures the electrical activity of the brain [START_REF] Zhang | Nonlinear electroencephalogram analysis of neural mass model[END_REF]. It has gradually gained considerable attention of neuroscientists due to its high spatial resolution. Also, functional Magnetic Resonance Imaging (fMRI) [START_REF] Kim | Cortical or Subcortical Neural Networks During Dynamic Neuromuscular Core Stabilization: A fMRI Blood Oxygen-Level Dependent (BOLD) Analysis[END_REF] is a useful technique for measuring different attributes of the brain activity. Indeed, both of EEG and fMRI techniques have shown their efficiency in revealing a better insight of the brain activity.

Despite that, both EEG and fMRI modalities suffer from some limitations. On the one hand, fMRI data analysis relies on slow mechanisms such as blood flow, oxygen consumption and metabolism that yield only an indirect measure of the evoked BOLD (Blood-Oxygen-Level Dependent) signals. Given this, it is difficult to estimate the neural activity only from hemodynamic response. Also, it has been clearly observed that fMRI suffers from an ill-posed temporal inverse problem (i.e poor temporal resolution). On the other hand, the EEG modality suffers from spatial inverse problem besides noise raised due to the unavoidable artifacts during the course of the experiment that can strongly degrade the EEG signals quality [START_REF] Mahjoory | Consistency of EEG source localization and connectivity estimates[END_REF] .

Both EEG and fMRI data measure neural activity. The EEG, revealing the event-related potentials (ERP) for investigating the psychophysiological states and information processing, is typical for the electrical activity. It is worth noticing that the evoked electromagnetic and metabolic responses detected by both EEG and fMRI are not essentially the result of the same primary neuronal processes [START_REF] Pflieger | Nonlinear analysis of multimodal dynamic brain imaging data[END_REF].

Given that strengths and weaknesses of fMRI and EEG complete each other, simultaneous EEG-fMRI acquisition applications are fast becoming a key technique to provide a more comprehensive understanding of the nature of the cerebral activity with its enhanced spatiotemporal resolution [START_REF] Ives | Monitoring the patients EEG during echo planar MRI[END_REF]. Thus, several works have been recently developed involving the EEG/fMRI fusion.

Mainly, the existing potential approaches devoted to develop EEG-fMRI integration procedures can be divided into two categories. In asymmetric methods, prior information of one modality is considered as a cause or a predictor of the second modality. In other word, it serves as a guide for the analysis of the other modality. Asymmetric approaches include: (i) fMRI informed EEG , known as fMRI-constrained EEG imaging, that employs spatial information from fMRI signal for reconstructing the source of the EEG signal [START_REF] Lei | A parallel framework for simultaneous EEG-fMRI analysis methodology and simulation[END_REF] and (ii) EEG-informed fMRI technique that is concerned with modelling the fMRI signal with features from EEG convolved with the Hemodynamic Response Function (HRF) [START_REF] Montes | Concurrent EEG-fMRI analysis by multiway partial least squares[END_REF][START_REF] Eichele | Prediction of human errors by maladaptive changes in event-related brain networks[END_REF]. EEG-fMRI symmetric methods are the second category of EEG-fMRI integration techniques. Mainly, the focus of these approaches is devoted to solving the inverse problem. They imply the use of a common forward or generative model to explore EEG and fMRI data jointly. Besides, they required some given knowledge from experience on the subject. [START_REF] Benar | EEG-fMRI of epileptic spikes concordance with EEG source localization and intracranial EEG[END_REF][START_REF] Daunizeau | Symmetrical event related EEG-fMRI information fusion in a variational Bayesian framework[END_REF][START_REF] Valdes Sosa | Model driven EEG-fMRI fusion of brain oscillations[END_REF][START_REF] Lei | EEG-fMRI fusion based on independent component analysis integration of data-driven and model-driven methods[END_REF][START_REF] Munck | The hemodynamic response of the alpha rhythm an EEG-fMRI study[END_REF].

To fuse the EEG and fMRI sources, a framework proposed in [START_REF] Daunizeau | Symmetrical event related EEG-fMRI information fusion in a variational Bayesian framework[END_REF] employs Bayesian theory to generate a common generative model. In this method, common parameters have been used to represent the spatial characteristics of both EEG and fMRI data. Therefore, the lack of detailed information about neurovascular coupling has not confused this fusion approach. Independent Component Analysis (ICA) has provided a powerful framework for integrating both modalities [START_REF] Lei | A parallel framework for simultaneous EEG-fMRI analysis methodology and simulation[END_REF][START_REF] Eichele | Mining EEG-fMRI using independent component analysis[END_REF]. To identify the component in EEG and fMRI independently, spatial ICA (sICA) for fMRI data and temporal ICA (tICA) for the EEG signal have been used. Due to the variation in spatial and temporal modalities, two methods have emerged: JointICA [START_REF] Calhoun | Method for multimodal analysis of independent source differences in schizophrenia combining gray matter structural and auditory oddball functional data[END_REF][START_REF] Calhoun | A review of group ICA for fMRI data and ICA for joint inference of imaging genetic and ERP data[END_REF] and parallel-ICA algorthim [START_REF] Eichele | Unmixing concurrent EEG-fMRI with parallel independent component analysis[END_REF]. Another similar approach is a parallel framework for Spatial-Temporal EEG/fMRI Fusion (STEFF) [START_REF] Lei | A parallel framework for simultaneous EEG-fMRI analysis methodology and simulation[END_REF]. Typically, STEFF has used an Empirical Bayesian (EB) model to combine temporal components from EEG and spatial components from fMRI obtained by ICA technique. Besides providing a mechanism for adjusting the problem of ignorance in ICA, STEFF has shown its ability to yield a flexible tool for imaging the EEG sources using information from fMRI data. For further detail, the reader is invited to see [START_REF] Mijovic | The why and how of JointICA results from a visual detection task[END_REF][START_REF] Mijovic | The dynamics of contour integration a simultaneous EEG-fMRI study[END_REF][START_REF] Poznanski | fMRI models of dendritic and astrocytic networks[END_REF]). Also, Dynamic Causal Model (DCM) is another framework for models of bioelectric and metabolic activity in neuronal populations [START_REF] Kiebel | Dynamic causal modelling of evoked responses the role of intrinsic connections[END_REF][START_REF] Chen | A dynamic causal model for evoked and induced responses[END_REF].

It is worth noticing that any symmetric fusion confronts three main problems: (i) The complexity of real metabolic-hemodynamic that renders the estimation of model-driven fusion problematic [START_REF] Valdes Sosa | Model driven EEG-fMRI fusion of brain oscillations[END_REF], (ii) The common substrate underlying EEG and fMRI signals should be identified by the fusion method and (iii) the uncertainty must be reduced by the fusion method when inferring this common subspace [START_REF] Daunizeau | Symmetrical event related EEG-fMRI information fusion in a variational Bayesian framework[END_REF]. Given this, an extensive research has been conducted on EEG-fMRI fusion data analysis. However, EEG-fMRI fusion requires a common theoretical framework to analyze and combine information sources with different natures (modalities).

To enhance the accuracy of the combined EEG-fMRI analysis, the current work introduces a new approach based on symmetric fusion using Dempster Shafer evidence theory (DS). Basic Belief function and combination rule have been employed as major keys in the proposed method. Mainly, DS theory permits to study different data like EEG and fMRI because the combination rule measures enable to join these data. BOLD signal is usually considered to have the same time evolution as the EEG [START_REF] Lei | EEG-fMRI fusion based on independent component analysis integration of data-driven and model-driven methods[END_REF]. Given this, the HRF model permits analyzing fMRI and EEG signals in simultaneous EEG-fMRI data. Thus, in this approach, the HRF model has been used as a key element [START_REF] Gotman | Epileptic networks studied with EEG-fMRI[END_REF][START_REF] Sturzbecher | Non-extensive entropy and the extraction of BOLD spatial information in eventrelated functional MRI[END_REF][START_REF] Sato | From EEG to BOLD brain mapping and estimating transfer functions in simultaneous EEG-fMRI acquisitions[END_REF][START_REF] Friston | Analysis of functional MRI time series[END_REF] for modelling both EEG and fMRI modalities by DS theory to fuse spatiotemporal information. To detect activated areas in the brain, the Potential-Based Hierarchical Agglomerative (PHA) [START_REF] Lu | PHA a fast potential-based hierarchical agglomerative clustering method[END_REF] have been employed. Therefore, PHA clustering method has time complexity of O(N 2 ) and allows to choose the max number of clusters , denoted as (C). PHA generates automatically a number of clusters that is less than or equal to the selected max number (C). Among these clusters, we have shown the target cluster and more information like the number of voxels by cluster and the centroid of cluster denoted by (c).

The remaining paper is organized as follows: Section two describes the basics of Dempster Shafer theory. Besides the brief review of the hemodynamic response model, the proposed model by using Dempster Shafer is introduced in section three. Section four gives details on the proposed approach. The conducted experiments on real and simulation data are discussed with the various findings in section five. At last, a conclusion is given.

2 Background of the Dempster-Shafer evidence theory Dempster-Shafer evidence theory has shown its efficiency in solving real problems [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF]. It has been widely used in medical imaging and signal processing like prediction, classification, segmentation and information fusion [START_REF] Liu | A new method based on Dempster-Shafer theory and fuzzy c-means for brain MRI segmentation[END_REF][START_REF] Cui | Dempster-Shafer multifeature fusion for pedestrian detection[END_REF]. In this work, two basic rules of DS theory are required: (i) Basic Belief Assignment (BBA) [START_REF] Shafer | A Mathematical Theory of Evidence[END_REF][START_REF] Yager | Dempster shafer belief structures with interval valued focal weights[END_REF], (ii) Belief function, bel (.) and (iii) Dempster rule of combination. Thus, these fundamentals are described respectively as follows:

Definition1: Basic Belief Assignment (BBA) is a key point of evidence theory, referred to as m(.), and has the following form:

m ∅ = 0, and m Ai !"⊆⊖ = 1 (1)
Definition 2: Belief function, bel (.).Let θ = θ ! , θ ! , … . , θ ! be a finite set of possible hypothesis. This set refers to the frame of discernment, and its power set is denoted by 2 where:

2 ! = ∅, θ ! , θ ! , … . , θ ! , θ ! ∪ θ ! , θ ! ∪ θ ! , … . θ bel (.) is described as: bel Ai = m Ai !"⊆! (2)
DS theory provides the combination rule that join several masses functions in one single masse probability of all modalities. The focus of this work has been also placed on employing the combined rule for pooling of evidence from two belief functions, Bel1 and Bel2

We assume two sources of evidence: B for EEG and C for fMRI are represented by respective subsets of θ : B1,...,Bm and C1,...,Cn. Probability density functions m1 and m2 are defined on B and C respectively. The Dempster rule of combination for combining two sets of masses m1 and m2 is defined as: m !" ∅ = 0 . A new value m12(A) is calculated for each subset of θ as follows:

m !" (A) = ! !!! m ! B m ! !∩!!!!∅ (C) (3)
where k denotes a measure of the amount of conflict between two evidences. If k = 1 the two evidences cannot be combined because their cores are disjoint. This rule is commutative, associative, but not idempotent or continuous. k is calculated as follows:

k = m ! B m ! !∩!!∅ (C) (4)
In what follows, the notations mfMRI (.), mEEG (.) and mEEG, fMRI (.) have been used to refer to BBA (fMRI), BBA (EEG) and the rule of combination respectively.

Modelling HRF by Dempster Shafer method

The canonical model proposed in [START_REF] Dempster | Upper and lower probabilities induced by a multivalued mapping[END_REF] is one of the numerous HRF models that have been developed in literature. It has a vital role to play in characterizing the onset of the stimulus. This work focuses on the study of the so-called sum two gamma functions [START_REF] Aguirre | The variability of human BOLD hemodynamic responses[END_REF]. The first is for modelling the peak whereas the second is for modelling undershoots. The HRF model that is typical for BOLD impulse response is presented in Figure 1.

Canonical hemodynamic response function model

To extract more information about the activated areas, the HRF function has been modelled as follows: First, HRF function has been divided into two hypotheses 𝜃 ! , 𝜃 ! . 𝜃 ! corresponds to modelling the neural activation and determining a peak (on activation). 𝜃 ! is assigned to modelling the undershoot (off activation). In particular, the first hypothesis is the focus of the current work. It is divided into two parts A and D where A stands for degrees of belief included in D. In other words, (D minus A) stands for uncertainty. Figure 2 depicts the proposed model. To compute basic belief assignments BBA and other Dempster Shafer metrics, the m(ti ) has been first computed using equations 5 and 6. Then, belief measures of each modality (fMRI and EEG) have been selected. Thus, the formula described bellow has been used, it consists in the transformation of each signal S(t) of fMRI data or EEG signal into density of probability function by using the formula:

m ti = ! ! S t dt !" !"!! for i = 1 to n m(ti)=0 , i=0 (5)
In the above equation, α is calculated as:

α = S t dt !" ! (6)
We obtain vectors of masses function for each modality. Then, these vectors have been employed for computing combination rule measures using the equations ( 3) and ( 4) that need to determine the intersection intervals 𝐴 = 𝐵 ∩ 𝐶 which are different of empty set ∅.

Proposed Method

As will be illustrated in Figure 4, the proposed scheme for EEG-fMRI symmetrical fusion is as follows: (i) pre-processing both fMRI and EEG data, (ii) generating HRF model by DS method, (iii) Calculating the mfMRI() for each voxel, (iv) computing mEEG() after extracting ERP from EEG signals and then, (v) mfMRI() and mEEG() metrics have been used to compute mfMRI, EEG() measure. This combination rule measure (mfMRI, EEG()) has been used to extract activated voxels via clustering. To separate the activated voxels from non-activated voxels, the clustering methods have been used because they can separate data into groups based on certain similarities. In this step, the PHA has been employed. As mentioned above, PHA generates automatically a number of clusters which is less than or equal to the selected max number (C). Furthermore, among these clusters, we show the target cluster and more information like the number of voxels by cluster and the centroid of cluster denoted by (c) as well. The process is described in Dempster Shafer algorithm (algorithm DS) and consists of six steps as follows:

Algorithm DS:

• IMPUT: fMRI data, EEG signals, HRF

• OUTPUT: m fMRI (.) , m EEG (.), m EEG,fMRI (.), activated areas

• step 1: first model HRF by using DS method, described in section 2

• step 2: project HRF model with fMRI signal and compute m fMRI () vector for each voxels using equations ( 5) and ( 6).

• step 3: from EEG to ERP

• step 4: the same process applied in step 2 has been repeated on ERP to compute vector of m EEG ()

• step 5: compute the combination rule measure using the equations aforementioned in ( 3) and (4)

• step 6: use PHA method to cluster vector of m EEG,fMRI (.) combination measure in order to identify active areas .

Experiments and results

The proposed method that relies on Dempster Shafer theory has been evaluated on both synthetic and real data. The performance of the introduced method has been based on statistical measurements. The MATLAB has been used for all implementations. At last, a comparison with jointICA has been provided to evaluate the efficiency of the proposed method. The evaluation metrics and experimental results are discussed below.

Evaluation metrics

This subsection describes the evaluating metrics serving to analyze the performance of the introduced approach and to establish a comparison with the jointICA. Basically, two main metrics have been given in [START_REF] Francois G Meyer | Spatiotemporal clustering of fMRI time series in the spectral domain[END_REF]: (i) True activation rate (TAR) presents the ratio between the number of voxels correctly identified as activated and the total number of truly activated voxels and (ii) false activation rate (FAR) is known as the ratio between the number of voxels incorrectly identified as activated and the total number of truly non-activated voxels. A more comprehensive description of these two metrics allows to correctly plotting the Receiver Operating Characteristics Curve ( ROC) [START_REF] Sorenson | ROC methods for evaluation of fMRI techniques[END_REF] for comparing between methods.

To compute these measures, it should be noted that jointICA has been performed on the basis of the t-test that uses a p-value or Z-threshold to map the activated voxels. The introduced method relies on the combination rule mEEG, fMRI ( ) to separate between activated and non-activated voxels. 𝑚 !"#$ (. ) 𝑚 !!" (. )

𝑚 !!",!"#$ (. )
Extracting ERP Extract active voxels by PHA clustering

Artificial dataset

First, the proposed method has been tested on artificial dataset developed by Correa et al. [START_REF] Correa | Comparison of blind source separation algorithms for fMRI using a new Matlab toolbox GIFT[END_REF] with their simulation algorithm. This dataset is available at (http://mlsp.umbc.edu/simulated_fmri_data.html). It contains eight sources where each simulated source is of 60 × 60 image with 100 point courses. However, in the current experiment, the virtual source for fMRI has been created into a two dimensional spatial map, as depicted in Figure 5, as the same size of the real fMRI data. Given that the information obtained from the brain is a stochastic process, a collection of sinusoid signals of the ERP data has been considered.

Let φ be the Gaussian random delay distributed with zero mean and unit variance and let d be the amplitude. Accordingly, ERP sources have been simulated by using the formula:

x n = 100×d ! sin 2πfn + φ , n = 1 … … 200 (7)
This sinusoid signal is with a mixture of very low frequencies from [0.014-0.811 Hz] The source S with a size of (8×3339) has been obtained by mixing the different FMRI sources. In the next step, the Demspter method has been applied on these data to get the mfMRI( ). Then, the same algorithm has been applied on ERP dataset to compute mEEG( ). Once both mEEG and mfMRI have been obtained, the mEEG. fMRI ( ) have been computed by using the combination rule. The obtained mEEG, fMRI ( ) has been clustered via PHA method. In this process, (k=8) has been chosen according to the number of simulated data sources. The obtained results are presented in Figure 5 where different slices have been clearly shown (mixed areas) in different sources. This figure presents the mixed activated areas detected by the proposed method. As it is clear in the first four slices on the left, the activated regions are depicted. The different sources (S1, S2,S5 , S6 and part of S4) figure 6 (a) . Slices on the right show the activated areas of each source in (S4, S6 and S8). Comparing our results with those obtained with jointICA, it can be seen that jointICA detects (S3,S4,S5,S6) sources (see Figure 6 (b)). Thus, the proposed approach has obviously shown its ability to yield a clear distinction between areas in terms of activation. Table1 describes the quantitative comparison between introduced method and jointICA using the true activation rate (TAR) and false activation rate (FAR). It's clear from this table that each method achieves better result more than 90.25%. Also from this table the proposed method detect mor accratly than jointICA for the sources (S1,S2, S6 ,S7, S8). Then jointICA outperform the proposed method for sources (S3,S4,S5). Table 2 show the average of measures (TAR and FAR).

AVG( TAR) AVG( FAR)

Proposed method 94,93% 0.051 jointICA 93.87% 0.061

From the outcome, it is very motivating to identify more activated areas. Thus, the introduced method can outperform better than the jointICA method in detecting more accurate activation for the simulated dat.

Real data

Besides simulation data, the current method has been evaluated also on a real fMRI and EEG dataset to access its performance. This data was described in [START_REF] Calhoun | Method for multimodal analysis of independent source differences in schizophrenia combining gray matter structural and auditory oddball functional data[END_REF], and was collected while participants were performing an auditory oddball. Prior to analysis, data goes through a series of preprocessing steps to identify and remove artifacts and validate model assumptions as well. Therefore, the fMRI slices require to be spatially realigned at first. However, spatial smoothing may cause unforeseen changes to occur into data. Thus, for better performance, spatial smoothing has been avoided. Then, the mean value has been subtracted from each of the time series and the variance has been normalized to a unit. The previous steps have been realized via SPM tools (Statistical Parametric Mapping).then the proposed method was applied on voxel. In this case, the mfMRI() and mEEG() was computed by using equations 5 and 6. Then, the both modalities fMRI and EEG data was fuse by employing the combination rule in order to get avector of mEEG, fMRI() by using equations 3 and 4 . Finally this vector is used for extracting activated voxels by using PHA clustering method. In addition we proceed to change the max number of cluster (k) by increment then we calculate the centroid (c) with the number of voxels in the target cluster. Then, the activated regions were shown. The same process has been repeated till no new cluster will be generated. Table 3 presents the results of mEEG, fMRI() clustering by PHA method. The plots depicted in Figure 7 present the various numbers of voxels by max number of clusters. By varying the max number (k), false positive voxels may be obtained if k is small or some voxels may be neglected if (k) is great. Thus, a threshold should be selected for better distinction of activated voxels. Thus, the average of the centroids of the target clusters has been used. By computing the average of the centroids (c), we get 0.006. Figure 8. (a) Illustrates slices that indicate activated regions while figure 8. (b) shows the activated areas for both components (red and blue). Qualitatively it is clear that introduced method detect activated region like jointICA method applied on the auditory data. In the term of quantitative comparison of the results obtained with the proposed approach and those obtained with jointICA that has been applied to analyze simultaneous EEG-fMRI data that is illustrated in Figure 9. As stated above in the current study, the proposed method depends on DS method that consists of two main components: The HRF model obtained by DS and mEEG, fMRI () measure. For selecting activated voxels, PHA method must be used for clustering this measure.

To validate the efficiency of the introduced approach, a comparison between the proposed algorithm and jointICA has been performed. It should be noted that jointICA is performed on the basis of the t-test that uses a p-value or Z-threshold to map the activated voxels. Given that the definition of p-value is similar to that of the false activated position rate in the jointICA, the pvalue can be used as a metric for choosing thresholds. Thus, several plots have been generated for both true and false activation rates at different thresholds (mEEG, fMRI, p-value). Consequently, the experiment has been performed by changing thresholds for computing the metrics. Figure 9 depicts the variation of true and false activation rates according to mEEG, fMRI() using the proposed algorithm. It seems obvious that the number of true activated voxels tends to be lower in true activation rate. However, the pace of false activation rate remains stable in the interval [0.001-0.005] then it gradually decreases. shows that the true activation rate remains the same till the p-value reaches 0.05 when it starts to fall down gradually. Figure 10 (b) illustrates that false activation rate continues to decrease gradually when the p-value increases and sometimes it remains fixed as in] 0.025, 0.035[. We then specifically compare the obtained results with those gained by using jointICA. The ROC curves for both methods are shown in Figure 11. It is obviously observed that the ROC curve obtained by jointICA is much lower than that obtained by the proposed method. While joint-ICA provided few true positives with few false positives, the new method has detected more true positives without increasing the number of false position. However, the area under ROC curve for DS method is 0.877 and 0.733 for jointICA method. Therefore, this comparison reveals that the proposed algorithm provides a better performance in terms of identifying more true positives than jointICA.

Conclusion

This paper introduces a new framework to identify activation areas in simultaneous EEG-fMRI data. Basically, the approach relies on Dempster Shafer theory taking into account the application of basic belief assignment and rule of combination. PHA method has been also employed to clustering the obtained measures for separate activate voxels from non-activated ones. An extensive performance analysis has been provided with both artificial and real datasets. Compared to joint ICA, the obtained activation map confirmed the efficiency of the proposed method for EEG-fMRI fusion in yielding a clear distinction between activated and non-activated areas. Although the focus has been placed on the analysis of auditory data, the introduced approach can be extended to explore various aspects of brain activity and to detect brain illnesses such as epilepsy.

  Figure 1.Canonical hemodynamic response function model

Figure 2 .

 2 Figure 2. : the proposed model of Dempster Shafer for Canonical hemodynamic response function

Figure 3 .

 3 Figure 3. Figure shows how to detect the interval A

Figure 4 .

 4 Figure 4. flowchart of the proposed method. The chart is devided into three main parts. The blue part concerns preprocess and modelisation of fMRI data.The green part is about EEG data and contains three steps; preprocessing, extracting ERP and modelisation by DS. The red part denotes clustering extract activated voxels

Figure 5 .

 5 Figure 5. Artificial data taken from [36]. (a) Artificial fMRI data contains eight sources where each simulated source is of 60 × 60 image with 100 point courses (b) the selected ERP simulation generated by using( equa. 7)

Figure 6 .

 6 Figure 6. (a) The obtained results by the proposed method that show the mixed area in each slice (b) the JointICA results

Figure 7 .

 7 Figure 7. (a) Variation in the number of voxels by number max of clusters k.

Figure 8 .

 8 Figure 8. : (a) Activated areas in some slices by proposed method (b) Activated areas in slices by jointICA method, provided in both coulor red and blue for each compenent

Figure 9 .

 9 Figure 9. . Number of true positives (a) and false positives (b) using different thresholds of (mEEG,fMRI).

Figure 10 .Figure 11 .

 1011 Figure 10.The number of (a) false activation rate (b) by using several p-values

Figure 10 (

 10 Figure10 (a)shows that the true activation rate remains the same till the p-value reaches 0.05 when it starts to fall down gradually. Figure10 (b) illustrates that false activation rate continues to decrease gradually when the p-value increases and sometimes it remains fixed as in] 0.025, 0.035[. We then specifically compare the obtained results with those gained by using jointICA. The ROC curves for both methods are shown in Figure11. It is obviously observed that the ROC curve obtained by jointICA is much lower than that obtained by the proposed method. While joint-ICA provided few true positives with few false positives, the new method has detected more true positives without increasing the number of false position. However, the area under ROC curve for DS method is 0.877 and 0.733 for jointICA method. Therefore, this comparison reveals that the proposed algorithm provides a better performance in terms of identifying more true positives than jointICA.

Table 1

 1 present true and false activation rates by the proposed method and jointICA

	source		S1	S2	S3	S4	S5	S6	S7	S8
	Proposed	TAR%	97.25	93.50	98.25	90.75	98.50	96.75	92.75	91.75
	method	FAR	0.0275	0.0650	0.0175	0.0925	0.0150	0.0325	0.0725	0.0825
	jointICA	TAR%	93.50	91.75	98.50	93.00	98.75	94.25	91.25	90.25
		FAR	0.0650	0.0825	0.0150	0.0700	0.0125	0.0575	0.0875	0.0975

Table 2 :

 2 Average of TAR and FAR of the proposed method and jointICA to the simulated dataset.

Table 3

 3 result of mEEG, fMRI () clustering by PHA method

	Number max of	Number of voxels	Center c
	cluster (C)
	3	2235	0.0053
	4	2184	0.0051
	5	2184	0.0051
	6	2184	0.0051
	7	2184	0.0051
	8	2184	0.0051
	9	1378	0.0062
	10	822	0.0069
	11	783	0.0068
	12	783	0.0068
	13	783	0.0068
	14	783	0.0068
	15	783	0.0068


S A Emine

G Bouleux

H Haouba

E Marcon
email: eric.marcon@insa-lyon.fr

Emergency Department Overcrowding Detection by a Multifractal Analysis

Keywords: Fractals Systems, Decision Support Systems, Detection, Signal Analysis, Data Flow Analysis

The observation of data which counts the arrivals at the ED in scales of hours or days reveals that these arrivals are characterized by the phenomenon of burstiness. The burstiness is a phenomenon that appears in the majority of Emergencies due to a batches of arrivals coming to the ED in a small time interval over a wide scales. The prediction of congestion in ED caused by the batched arrivals seems important to the medical staff. The proper modeling data like arrivals process that has this scaling property is relevant with the self similar time series. The goal of our study is to establish that the time series of patient flow displays fractal behavior, whose quantitative characteristics vary with time. To discover whether a Patients flow has begun to congest, we analyze the time series of patients arrivals data with a Multi-Fractal Detrended Fluctuation (MF-DFA) algorithm.

INTRODUCTION

Overcrowding in the Emergency Department (ED) is a worldwide problem impairing the ability of hospitals to offer emergency care within a reasonable time frame [START_REF] Lin | Estimating the waiting time of multi-priority emergency patients with downstream blocking[END_REF]; [START_REF] Bouleux | Early index for detection of pediatric emergency department crowding[END_REF]. It is a reflection of several internal or external effects such as the capacity of the ED or the stream of arrivals. During the last several years, many research papers have been devoted to the topic of patient flow or stream of arrivals in the ED. One can see for example, the following manuscript and the references that are quoted in it [START_REF] Armony | On patient flow in hospitals: A data-based queueing-science perspective. an extended version[END_REF]. The work proposed in [START_REF] Monte | On the selfsimilar distribution of the emergency ward arrivals time series[END_REF] analyze the ED admissions time series and search the distribution that best models it. In her thesis, [START_REF] Jlassi | Amélioration de la performance par la modélisation des flux logistiques des patients dans un service d'urgence hospitalier[END_REF] improved the performance of patient logistics flows in a hospital emergency department by applying to the hospital the analysis methods and the resolution tools from the manufacturing field. By using the time series, specifically an ARIMA model, [START_REF] Kadri | Modélisation et prévision des flux quotidiens des patients aux urgences hospitalières en utilisant l'analyse de séries chronologiques[END_REF] modeled the daily flows of patients for the pediatric emergency of Lille, France, Hospital. Finally, the most commonly encountered model for the hospital emergency arrivals is definitely the Poisson process. This modeling is for example introduced in [START_REF] Lin | Estimating the waiting time of multi-priority emergency patients with downstream blocking[END_REF] where the authors have estimated the waiting time of multipriority emergency patients using the Fast track system and a queueing model. They modeled the arrival process of demand accessing the ED by a Poisson process so. In [START_REF] Koizumi | Modeling patient flows using a queuing network with blocking[END_REF], the authors divide the whole day into several segments such that the arrival rate in each segment can be assumed to be constant and affirme that the arrival process of emergency patients in most time segments can be modeled as a Poisson process with varying rates. Nevertheless, these approaches turn out not to work for ED; the pattern of arrivals found by the staff of the emergency ward service suggested that this model was inadequate [START_REF] Monte | On the selfsimilar distribution of the emergency ward arrivals time series[END_REF]. Indeed, the time series of arrivals in ED exhibit burstiness over a wide range of time scales due to pikes of arrivals process batched over hours, days or weeks. A process which takes into account features related to bursts is known to be a self-similar process with the Long Range Dependence (LRD) property [START_REF] Radev | Modeling and simulation of self-similar traffic in wireless ip networks[END_REF]. Such a process has a statistical structure that repeats itself on subintervals or time scales [START_REF] Espen | Introduction to multifractal detrended fluctuation analysis in matlab[END_REF]. This statistical characteristic has shown its good performances in network traffic congestion prediction [START_REF] Kim | Congestion prediction of self-similar network through parameter estimation[END_REF] and some works initiated the reflexion on ED admissions modeling by self-similar processes [START_REF] Monte | On the selfsimilar distribution of the emergency ward arrivals time series[END_REF]. If the self-similarity propose a global scaling description of the process, the fractality and multifractily are properties which focus on local scaling behavior and appears as key features for congestion detection as well [START_REF] Feng | A better understanding of long-range temporal dependence of traffic flow time series[END_REF]. Consequently, we propose in this paper a Multifractal analysis of the ED arrivals time series with the aim of being able to detect the burstiness and congestion caused by this phenomenon. We present next the basic notions of our study and recall the fundamental of the self-similarity properties for a time series. Section 3 is dedicated to the local study of time series exhibiting scaling by the multifractal analysis. We illustrate then this analysis on real ED data admissions flow, in view to depict the congestion caused by the burst Preprints of the 10th IFAC Symposium on Biological and Medical Systems São Paulo, Brazil, September 3-5, 2018 Mo6Track B.2 10th IFAC Symposium -BMS 2018, São Paulo, Brazil, September 3-5, 2018 of arrivals at the ED. Finally, a conclusion is drawn in the last section.

FRACTAL ANALYSIS

Notion Of Self-Similarity

A time series X t , t ∈ N with mean µ(t), variance σ 2 (t) and autocorrelation function ρ

(k), k ≥ 0 is called a stationary process if E[X t ] = µ(t) = µ ∀ t ∈ N (1) with E[.] the mathematical expectation, and cov(X t , X t-j ) = E[(X t -µ)(X t-j -µ)] = γ j , ∀ t, ∀ j. (2)
For each m = 1, 2, . . . , constructs from the covariance stationary process X t a new covariance stationary time series X (m) = {X (m) k : k = 1, 2, ...} with corresponding autocorrelation function ρ (m) , given by

X (m) := 1 m (X km-m+1 + ... + X km ), k ≥ 1. ( 3
)
X t is said to exhibit self-similar (or asymptotically self similar) property, if for all t = 0, 1, 2, ...,

X t d = m (1-H) X (m) (4)
holds ∀m (or as m → ∞), where 1 2 < H < 1 and d = is the finite dimensional law equality. X t is said to exhibit (exactly) second-order self-similar property with self-similarity parameter

H if m (1-H) X (m) t
has the same autocorrelation function as X t , ∀ m and ∀ t.

cov(X (m) ) = σ 2 m 2(H-1) and ρ (m) (k) = ρ(k), k ≥ 0. (5)
In other words, X t is exactly second-order self-similar if the aggregated processes X(m) are indistinguishable from X t with respect to their first and second order properties Strza lka et al. ( 2012). In addition, X t is said to exhibit asymptotically second order self-similar with selfsimilarity parameter

H (H-a.s.o.s.s.) if ∀ k, lim m-→∞ ρ (m) (k) = ρ(k) (6)

Notion of Long Range Dependence

The Long Range Dependence (LRD) as a phenomenon was first observed by [START_REF] Hurst | Methods of using long-term storage in reservoirs[END_REF] when he studied the flow of water in the Nile river. Many works have emerged since but [START_REF] Mandelbrot | Une classe de processus stochastiques homothetiques a soi-application a la loi climatologique de hurst[END_REF] is certainly the reference papers. Mathematically, a stationary process X t is said to exhibit short range dependence if its autocorrelation function ρ(k) satisfies

∞ k=1 ρ(k) < ∞. (7)
Conversely, if the autocorrelation function ρ(k) is not summable, the process X t is said to exhibit Long Range Dependence such as

k≥0 |ρ(k)| = ∞. (8)
An equivalent definition is often proposed regarding whether the autocorrelation behaves like a power function decaying to zero hyperbolically

ρ(k) ∼ L ρ (k)k -α 0 < α < 1 (9)
with L ρ (k) a slowly varying function [START_REF] Bingham | Regular variation[END_REF]. Lastly, the Karmata definition of this class of functions is the set of functions such that the following equation holds lim

x-→∞ L ρ (k)(ax) L ρ (k)(x) = 1. ( 10
)

The connection of Self similarity and LRD

The connection between Self-similarity and long memory phenomena was studied extensively by [START_REF] Samorodnitsky | Long memory and selfsimilar processes[END_REF]. It was shown that asymptotically second order selfsimilarity property implies LRD propertyLazarou et al. ( 2009); [START_REF] Park | Self-similar network traffic and performance evaluation[END_REF]. The autocorrelation of a Self Similar process with Hurst exponent H and finite second order momentEmbrechts and Maejima ( 2000) is then given by

ρ(k) = 1 2 ((k + 1) 2H + (k -1) 2H -2k 2H ). (11)
From this identity we straightforwardly obtain

ρ(k) ∞ ∼ H(2H -1)k -2(1-H) (12
) in asymptotic sense so. If 0.5 < H < 1, then the process admits LRD, we say it is persistent. On contrary, when 0 < H < 0.5 it is short range dependence and is qualified as antipersistent [START_REF] Park | Self-similar network traffic and performance evaluation[END_REF]. The Hurst exponent H is also understood in that context as a measure of the decay speed of the autocorrelation function tail. It captures then both, the intensity of LRD Vemuri (2005) and the degree of self-similarity Cervantes-De la [START_REF] Cervantes-De La Torre | Fractal dimension algorithms and their application to time series associated with natural phenomena[END_REF].

The Fractal Analysis and the ED arrivals process

The burstiness phenomenon we have mentioned in the introduction is very much related to the Hurst exponent. To describe quantitatively the burstiness due to the batches of arrivals processes, which depend on the time scale, Gusella [START_REF] Gusella | Characterizing the variability of arrival processes with indexes of dispersion[END_REF] propose the Index of Dispersion for Counts (IDC). It is defined as

IDC(N ) := var( N t=1 X t ) E( N t=1 X t ) (13)
where X t is the self-similar arrivals process at the ED. In the case of persistence for X t , the IDC varies exponentially with the interval length and follows the relation

IDC(N ) ∼ c N 2H-1 (14) log(IDC) ∼ (2H -1) log(N ) + const. ( 15
)
10th IFAC Symposium -BMS 2018, São Paulo, Brazil, September 3-5, 2018 characteristics vary with time. Such behavior is called multifractal.

THE MULTIFRACTAL ANALYSIS

The multifractal framework has become a major concept and tool involved in practical analysis of scaling in data.

The multifractal (MF) framework relates the scaling properties to the local regularity properties [START_REF] Lashermes | Comprehensive multifractal analysis of turbulent velocity using the wavelet leaders[END_REF]. In general, multifractality in time series is due to different LRD of the small and large fluctuations in the time seriesKantelhardt et al. (2002).

A self-similar stochastic process {X t : t ∈ [0, T ]} is multifractal [START_REF] Thompson | Multifractal detrended fluctuation analysis: Practical applications to financial time series[END_REF] if it satisfies

E[|X t |] q = c(q)t τ (q)+1 , ∀ t ∈ G and q ∈ Q (16)
where 0 < T < ∞, G and Q are open intervals on [0, T ] ⊂ G and [0, 1] ⊂ Q. The function τ (q) is called the scaling function [START_REF] Calvet | Multifractality in asset returns: theory and evidence[END_REF] of the multifractal process. The Multifractal Formalism is based on the calculation of two sets of coefficients: the Holder exponents h which quantifies the local regularity of a process X t and the MultiFractal Spectrum that quantifies the multifractality of the process [START_REF] Xiong | The timesingularity multifractal spectrum distribution[END_REF].

The Holder Exponent and the Spectrum of Singularities

The local regularity of the sample path of stochastic process X t is commonly studied via the notion of pointwise Holder exponent. Recall that X t is assumed to be locally bounded, i.e it belongs to C α (t 0 ), t 0 ∈ R with α ≥ 0, if there exists a constant C > 0 and a polynomial P satisfying deg(P ) < α and such that, in a neighborhood of

t 0 |X t -P (t -t 0 )| ≤ C|t -t 0 | α . (17)
The Holder exponent of X t at t 0 is therefore defined as

h(t 0 ) = sup{α : X ∈ C α (t 0 )}. ( 18
)
The value of the Holder exponent is interpreted as follows: the closer to 0 h(t 0 ) is, the more irregular (or singular) at point t 0 the function is. In contrast, a larger value of h(t 0 ) is related to a smoother (more regular) behavior at t 0 [START_REF] Lashermes | Comprehensive multifractal analysis of turbulent velocity using the wavelet leaders[END_REF]. The points t 0 associated with a specific Holder exponent value h are distributed on interwoven fractal subsets E(h) such as

E(h) = {α : h(t 0 ) = α}. (19)
The fluctuations, along time t of the Holder exponent h are usually described through the singularity (or multifractal) spectrum, labelled D(h) and defined as the Hausdorff dimension [START_REF] Przytycki | On the hausdorff dimension of some fractal sets[END_REF] of the set of points E(h) where the Holder exponent takes the value h, this leads then to

D(h) = dim H (E(h)). (20)
The MultiFractal analysis aims at characterizing or derive the spatial distribution of Holder exponents with singularity spectrum D(h) [START_REF] Lashermes | Comprehensive multifractal analysis of turbulent velocity using the wavelet leaders[END_REF]. Some relations between the two sets of multifractal scaling exponent have been established. τ (q) = qh(q) -1 (21) with h(q) the generalized Hurst exponent or the q-order Hurst exponent (F q (s) ∼ s h(q) ) and τ (q) the scaling exponent and (Z q ∼ s τ (q) ) the generalized multifractal dimensions. We obtained consequently the following relation

D(q) = τ (q) q -1 = qh(q) -1 q -1 . (22)
Next, the spectrum of singularity is related to the scaling exponent via the so called Legendre transform α = τ (q) = h(q) + qh(q), (23) and

f (α) = qα -τ (q) = q[αh(q)] + 1 ( 24
)
where α is the Holder exponent, while f (α) denotes the dimension of the subset of the series that is characterized by α. When Holder exponent takes a unique value H at every point t 0 , the multifractal process becomes a monofractal one and we obtain

f (α) = 1 if α = H 0 if α = H and α > 0 (25)

The MF-DFA Method

We propose here to introduce the Multifractal Detrend Fluctuation Analysis, originaly introduced by [START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF] as a generalization of the DFA method. It has a practical interest in giving singularities information, reason why we give its implementation principles in what follows. Suppose that X t , t = 1, . . . T is a series of the arrivals by days with mean µ. Timely integrate the detrended series X t to obtain

Y t = t k=1 (X k -µ) (26)
Next, divide the profile Y t into N s = N/s (with . the integer part ) non-overlapping segments of equal length s. Then, calculate the local trend for each of the 2N s segments by a least-square fit of the series and determine the variance with

F 2 (ν, s) := 1 s s i=1 {Y [(ν -1)s + i] -y ν (i)} 2 , ν = 1, ..., N s and F 2 (s, ν) := 1 s s i=1 Y [N -(ν -N s )s + i] -y ν (i) 2 , ν = N s + 1, ..., 2N s
where y ν (i) is the fitting polynomial in segment ν. Compute the average over all segments ν to obtain the q-th order fluctuation function The q-order singularity exponent

The q-oder singularity dimension -40 -20 0 20 40 0.7 0.8 0.9 1.0 q-order The q-order Hurst exponent

Figure 5. The q-order Hurst exponent of Arrival in ED estimated by the MF-DFA.

F q (s) = 1 2N s 2Ns ν=1 (F 2 (s, ν)) q 2 1 q . ( 27
)
If the time series X t has a long-range power-law correlated, then F q (s) increases for large values of s as a power-law which means that F q (s) ∼ s h(q) .

The Arrivals Process and the MultiFractal Analysis

In this section we analyse the time series of ED arrivals of Saint Etienne University Hospital (CHU), France, from 06/01/2013 to 08/31/2017. To estimate the q-order Hurst exponent of Figure 5, the local Hurst exponent Figure 6 and the the MultiFractal Spectrum (figure 4), we use the MF-DFA method [START_REF] Kantelhardt | Multifractal detrended fluctuation analysis of nonstationary time series[END_REF]. The Mul-tiFractal Spectrum (MFS) of the time series is plotted in Figure 4. The MFS is a downward concave function with maximum value 1 and minimum value around 0.45. This particular shape justifies the multifractal approches.

To depict the congestion in the stream, we observe the local Holder exponent through Figure 6. A cursory glance at the Holder exponent plot which depicts the temporal dynamics, shows that the arrivals time series admits a high variability of corresponding so to fractal behavior. However, it is not possible to distinguish a regularly regime of tensions. The average of arrivals at ED is 136 patients a day. Overall, the local Hurst exponent oscillates between 0.6 and 0.8 and when it exceeds 0.8 we notice that the arrivals are globally below the daily average of arrivals. On the other hand, when it is down 0.6, the system registers peaks of arrivals more than average.

CONCLUSION

Led by recent statistical analysis taking into account burstiness effects in view to early predict congestion on computers networks and traffic networks, we have proposed in this work to apply the underlying fractal methods to the stream of patients' arrivals at an Emergency Department (ED). The stream of arrivals at ED is characterized by a high variability not only on the dynamics of the process but also on strong irregularities of its scaling exponent. To fully characterize the arrivals flow, the fractal methods appeared to be insufficient and MultiFractal Analysis (MFA) is thus recommended. MFA only based on patients' arrivals at the ED gives information on the regimes of the arrivals, but it does not reveal all the system characteristics, particularly the burstiness induced by capacity constraints. Thus, MFA could be giving more information if the Length Of Stay (LOS) of patients at the ED were considered and used in conjonction with the patients' arrivals.

Figure 1 .

 1 Figure 1. The cumulated time series of Arrivals, (a) time series aggregated over 1 hour (black), (b) time series aggregated over 6 hours (red) and 8 hours (blue).

Figure 2 .Figure 3 .

 23 Figure 2. Plot in logarithmic scales of the Index of Dispersion of Counts (IDC) for the ED arrivals process with respect to hours (red) and days (blue) scales.The Hurst exponent of both is H=0.54The IDC can also served as a statistical method for testing the self-similarity as illustrated by Figure2. This plot gives the picture of IDC for the ED arrivals for different aggregation sizes N . The Figure indicates that the hurst exponent for the ED arrivals process is around H ∼ 0.54.Another popular approach adopted for self-similarity characterization is the Detrended Fluctuation Analysis (DFA), originally introduced by[START_REF] Peng | Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series[END_REF] and has

Figure 4 .

 4 Figure 4. The MF spectrum of The time series arrivals in ED ever days estimated by the MF-DFA.

Figure 6 .

 6 Figure 6. The local Hurst exponent of the arrivals by day estimated by the MF-DFA.

10th IFAC Symposium -BMS 2018, São Paulo, Brazil, September 3-5,