Context and method

Gallica (http://gallica.bnf.fr) is one of the major digital libraries available for free via the Internet. It provides access to 4.3 million documents of various types: printed documents (books, press and magazines), manuscripts, sound and iconographic documents, maps and plans, etc. Gallica receives around 1.5 million visits per month.

In the context of the Bibli-Lab research partnership between the BnF and Télécom ParisTech, and with the support of TeraLab, a new analysis of Gallica servers' connection logs was carried out, applying machine-learning methods to them. The aim was not to collect information on users or their profiles but rather to use logs 1 , which act as records of usage, as a basis for identifying typical clickstreams. For 15 months (April 2016-July 2017), Adrien Nouvellet, a researcher on a postdoctoral contract and under the supervision of four of Télécom ParisTech's research professors 2 , developed a data clustering algorithm enabling grouping of Gallica sessions with similarities in sequencing and duration of actions 3 . Logs analysed covered a range of time frames, from a week to a month, with systematic checking of the stability of models obtained.

Such learning methods take advantage of the very factor that undermines traditional methods for gathering information on usage: the extremely high numbers of connections (45,000 visits per day on Gallica). Such numbers call the representativeness of online surveys into question -representative first and foremost of the most committed internauts, but not of all internauts.

Despite the power of the algorithms involved, machine learning also requires numerous decisions to be taken, necessitating availability of other sources of knowledge on usages and users. For this reason, the preferred methodological choice was to have statistical models interact with results obtained from other approaches 1 Files containing all requests received by servers. Among other information required for knowledge of usage, logs contain the IP address (single connection identifier, anonymised for the purposes of this project), date and time (to the nearest second) of the request, the user's provenance (referring website), and the http request, which, in cases of requests for a Gallica document, contains its permanent identifier ARK. 2 Florence d'Alché-Buc and François Roueff from the Image, Data and Signal (IDS) Department; Valérie Beaudouin and Christophe Prieur from the Economic and Social Sciences (SES) Department. 3 The five "actions" identified in the logs are consultation of the homepage, consultation of collections guides (presentation of collections and blog), use of the internal search engine, consultation of a document in the Gallica interface, and downloading.

 (ethnographic observations, interviews, etc. 4). Such dialogue enabled the researchers involved to: a) set departure parameters (definition of a session and the elementary actions composing it); b) check models obtained, which were highly sensitive to technical artefacts; and c) propose initial keys to interpretation.

The interest of the work carried out on the Gallica logs persuaded the BnF and Télécom ParisTech to add a further stage to the research, this time lasting four months (July-November 2017) and devoted to Data BnF logs well as clickstreams between Gallica, Data BnF and the BnF General Catalogue. Available online since July 2011, the Data BnF project (http://data.bnf.fr) endeavours to make the data produced by the BnF more useful on the Web. This website gathers BnF resources and external resources on HTML pages devoted to an author, a work, a subject, a year or a place. Using Semantic Web technologies, these pages are automatically created with data and identifiers that are in the different databases of the library.

The main results of the two components of this research are presented below. The full report of the Gallica log analysis, including a presentation of the methodology, is available in HAL 5 .

Part 1 / Gallica: main results

Large numbers of very short sessions among Gallica users

Typical of life on the web today, most Gallica users make very brief visits to the site: 50% take less than 12 seconds; 30% only make one request and only 8% of sessions consult more than four separate documents 6 . In addition, only one in every four sessions makes use of the internal search engine, a fact that has two possible explanations: a) a high percentage of visits are not made in order to explore collections (but rather for a single consultation of a specific document); and b) a number of Gallica users substitute Google for the internal search engine, not just at the beginning of but also during a session, when they want to start a fresh search.

The study of referring sites (address of the web page at the origin of a request) shows that Gallica users' web provenance has an influence on session "depth", measured here by the number of actions carried out on Gallica. Although, unsurprisingly enough, Google is the main referring site irrespective of how many actions take place during a session. The same is not true of Facebook: the social network is best represented at the origin of sessions comprising between 2 and 4 actions (30% of sessions). In contrast to received ideas then, sessions originating with Facebook are not necessarily simple "one-click" sessions. For sessions comprising over four actions (40% of sessions), referring sites -after Google but before Wikipedia and Facebook -include the bnf.fr domain and a first theme-based site focusing on genealogy (Geneanet). This is proof of the importance of theme-based sites (which attract a large proportion of amateur researchers 7 ), at the origin of "in-depth" consultations of Gallica.

"Atypical" clickstreams

In order to have sufficiently rounded models available, a method for forming "clusters" (groups of similar sessions) was applied to sessions comprising more than five actions -only 35% of all Gallica sessions.

Example 1: 15-category classification

Although website design always supposes its "normal" use (e.g. homepage > internal search engine > consultation of document), clusters testify to the wide variety of clickstream patterns in Gallica consultations. In the first model of clusters obtained, where only sequences of actions are taken into account, 53% of sessions correspond to sequences 4 See below for references to other studies carried out upstream of or in parallel to this research. 5 Nouvellet A., Beaudouin V., D 'Alché-Buc F., Prieur C., Roueff F. (2017), "Analyse des traces d'usage de Gallica : Une étude à partir des logs de connexions au site Gallica", Rapport de recherche, Télécom ParisTech, Bibliothèque nationale de France, online: <https://hal.archives-ouvertes.fr/hal-01709264>. 6 A short session should not necessarily be seen as a failure, however: it might be a case of an internaut finding what he/she is looking for immediately by clicking on the right link (hypertext). 79% of single-task sessions on Gallica consist of a simple document request 7 Beaudouin V., Pehlivan Z. ( 2016), "Cartographie de la Grande Guerre sur le Web" (Cartography of the Great War on the Web), Final Report on Phase 2 of "The Future of Online Digitized Heritage: The Example of the Great War" project, online: <https://hal.archives-ouvertes.fr/hal-01425600/document>.

of simple document consultations that did not make use of the homepage, did not download and did not use the search engine. Collection presentation pages appear in a single small (2.5%) cluster, proof that such pages are not on most Gallica users' routes; their consultation constitutes a form of behaviour distinct from all the others observed.

Incorporation of the time factor provided a more accurate idea of sessions as more than simple sequences of actions: it brings together sessions that were previously assigned to different clusters but which turned out to have the same temporal "outline". As an example, consultation time for a "view" of the same document on Gallica may vary from 0 to 52 minutes, depending on conventions chosen for analysis 8 . It would seem reasonable not to assimilate such heterogeneous levels of engagement.

With the new model, the largest cluster (42%) brings together sessions taking an average of 7 minutes and alternating consultation actions (3 minutes on average) and search engine actions (2 minutes on average), with more random and briefer presence of other actions -a model that is therefore closer to that of the site's designers. Other clusters also deserve attention as they enable quantification of behaviour patterns that are both simple and typical: 28% of sessions are dominated by downloading activity, often associated with consultation; 13% of sessions only comprise alternating consultation and search engine actions, with more time devoted to the former; in contrast, 3.5% comprise long search sequences (11 minutes) terminating with brief consultations. Finally, although use of the homepage during a session usually takes no longer than 30 seconds, 7% of sessions spend an average of 4 minutes on it (not necessarily continuously) within a varied sequence of actions that even occasionally includes consultation of collection presentation pages. This final result shows that, even though it may scarcely be apparent, there are people who consult the homepage and know how to make good use of it for careful and varied consultations. The largest clusters should not be allowed to conceal forms of behaviour that are non-negligible at a quantitative level and of interest as regards development of audiences and uses.

If we now return to all sessions, however many actions they include, it can be seen that 4% of them never get past Gallica's homepage (13% of single-action sessions) and do not actually enter the site itself.

Little diversity in types of documents consulted during a session

In order to measure the diversity of documents consulted, logs were provided with the metadata of documents contained in the Open Archives Initiative (OAI 9 ) repository. Correcting the assertion made in the 2016 online questionnaire on this point 10 , it is clear that daily-press publications are the type of document most often consulted, ahead of monographs and images, even though, compared to the number of press titles available in Gallica, the consultation ratio is the lowest. This is surprising: despite the exploration facilities provided by the web's interfaces and the much-vaunted role of serendipity, Gallica consultations are still largely single-type. This is the case with 45% of sessions in which more than 5 documents are consulted, with predominance, as might be expected, of sessions focusing on daily-press publications and monographs alone. Sessions analysed, although longer than average, bear witness to users consulting the site in "silo" fashion (as if organisation of collections and research practices were still compartmentalised), as was shown by the 2012 study of document requests in the BnF's "Rez-de-jardin" reading rooms 11 . One challenge for the Gallica interface will be to encourage "bouncing" from one type to another (from a manuscript by Apollinaire to a recording of his voice, for example). Only 3% of sessions covering more than 5 documents explore almost all document types.

The most frequent actions carried out during a session vary greatly depending on documents requested. During sessions exclusively devoted to press publications, Gallica users spend most of their time consulting (i.e. making requests containing an ARK) and much less of it downloading. Unlike sessions devoted to manuscripts, where downloading takes up most of users' time 12 . This type of analysis can help prioritise developments of functionalities by type of document. Tools for exploring press collections and designed to facilitate refinement (e.g. search by named entities) would therefore appear to be necessary. 8 The video ethnography carried out in parallel to this analysis evidenced the existence of very lengthy consultations of single views (see Rollet N., Beaudouin V., Garron I. ( 2017), "Vidéo-ethnographie des usages de Gallica" (Video ethnography of uses of Gallica), Final Report on Phase 2, online: <https://hal.archives-ouvertes.fr/hal-01709210>), which led to revision of the definition of a session: a session on Gallica ends when the time between two requests exceeds 60 minutes (where, for other web services, received practice advises 10 minutes). We should bear in mind, however, that logs tell us nothing about user activity outside Gallica. 9 The Open Archives Initiative (OAI) is an organization to develop and apply technical interoperability standards for archives to share catalog information (metadata). 10 TMO Régions (2017), "Enquête auprès des usagers de la bibliothèque numérique Gallica" (Survey of Gallica Digital Library Users), survey report, online: <http://www.bnf.fr/documents/mettre_en_ligne_patrimoine_enquete.pdf>. 11 Pardé Th. (2015), "Les usages documentaires dans une bibliothèque de Recherche" (Documentary usages in a Research Library), Bulletin des Bibliothèques de France (BBF), no.5, p.112-119, online: <http://bbf.enssib.fr/consulter/bbf-2015-05-0112-002>. 12 This result is clarified by qualitative analysis: interviewees stated that they consulted the press in Gallica as its reading tools were well suited to the purpose (zoom functions in particular), whereas they downloaded monographs, which could be read just as easily offline (see Rollet N., Beaudouin V., Garron I. ( 2016), "Je pars d'un sujet, je rebondis sur un autre: pratiques et usages des publics de Gallica" (I start off on one subject and bounce to another: Gallica users' practices and usages), exploratory qualitative study, Final Report on Phase 1, September 2016, online: <https://hal.archives-ouvertes.fr/hal-01709238>). The 2016 survey by TMO Regions also drew attention to the increase in careful consultation in Gallica users' stated practices (TMO Régions, 2017).

Social networks: the impact that methods of publication on Facebook have on direct users

Communication on social networks has a proven impact on Gallica users, and, as we have seen, Facebook is well represented among referring sites. Analyses confirmed this with regard to a specific event: the peak period for consulting Voltaire's works on the anniversary of his death in 2016 was compared with the success of a publication on Gallica's Facebook page the same day (1,055 likes and 568 shares). Furthermore, although numbers of publications per week on the page have remained relatively stable since its creation, numbers of reactions per publication have increased significantly.

A study on types of links to Gallica included in publications has shown that they affect the number of "clicks": an active link in an image generates 25 times as many visits to Gallica than an active link in the text (with indication of the URL). This result led the team responsible for the Facebook page to modify its publication methods.

Part 2 / Data BnF and the General Catalogue

The aim of the Data BnF project (<http://data.bnf.fr>) is to make the BnF's data more useful on the Web. It involves transforming existing data, enriching and interlinking the dataset with internal and external resources, and publishing HTML pages for browsing by users and search engines. The raw data is also available in RDF following the principles of linked open data architecture Data BnF consultations follow a particular pattern as far as times of day are concerned: peak times for consulting the interface occur during the day (at 11 a.m. and 3-4 p.m.), whereas the peak period for Gallica is during the evening (between 8 and 10 p.m.). In-depth visits to Data BnF may also be quantified: 800 sessions a day include over 4 actions (5% of the total) and 250 sessions (1.5%) consult more than 4 entity pages. "Auteur" is by far the most consulted entity, although "Thèmes", "Spectacles" and above all "OEuvres" (the second most consulted) are consulted proportionally more often than their presence in Data BnF would suppose.

Over the course of a month, 12% of two million authors were consulted -a considerable figure. The function providing a breakdown of authors consulted at least once in a month enables formulation of the hypothesis (which must nonetheless be tested out over a longer time period) that most authors will have been consulted after enough time has passed and that none of them really leads the others. There would therefore be no "dark zone" in Data BnF, whose listing of all webpages seems to function as it should.

As initially desired by the Metadata Department, clickstreams between Gallica, Data BnF and the BnF General Catalogue were then subjected to analysis with a view to understanding how today's internauts access BnF documents and metadata on the web: how they enter applications and travel -or otherwise -from one to another. The large number of consultations of Gallica alone results in only 4% of all sessions combined on the three sites being multisite. Among such multisite sessions, Data BnF is by far the most frequently in evidence (87%), ahead of Gallica (69%): clear proof that it is successful in playing its role as a "pivot" between the BnF's main documentary and bibliographic services.

Data BnF acts as a gateway to Gallica documents before acting as one to data contained in the BnF General Catalogue. Half of all bounces are between Data BnF and Gallica, with only one out of every three bounces between Data BnF and the catalogue. Direct access to documents is therefore the prime motive among users who, after visiting Data BnF, continue their search in the bnf.fr domain. Such bounces sometimes occur unexpectedly from Gallica to Data BnF or, more rarely, from the BnF General Catalogue -probably due to the browser's navigation tabs, as there are no direct links from Gallica documents or the BnF General Catalogue to Data BnF.

Grouping multisite sessions into "clusters" confirms that:

• Data BnF is not simply a gateway to other services: in over a quarter of multisite sessions, lengthy sessions on Data BnF precede resulting consultations of Gallica (13%) or the BnF General Catalogue (5%);

•

10% of multisite sessions take full advantage of all three services and may, by general agreement, be referred to as "expert";

•

The most "expected" sorts of sessions -brief consultations of Data BnF that switch to long consultations of the BnF General Catalogue (initial idea of a "documentary pivot") -actually account for only 3.5% of multisite sessions.

Finally, as with Gallica, most Data BnF sessions (54%) only consult a single entity type, most often "Auteurs". Bitype sessions are usually sessions in which one type is consulted for a relatively long time before a brief switch to the second, but with no real alternation of different types during a session ("Auteurs" / "OEuvres" alternations only occur in 1% of sessions).

Avenues for work

Measuring evolutions in usage as the interface evolves

For Gallica, all processing was carried out on logs prior to 15 June 2016, i.e. before the redesign of the homepage. Renewing such analysis would enable very direct measurement of the impact that the evolution of the interface has on behaviour patterns.

Exploring segments of the public important to development of the interface

Analysis of logs enables identification and quantification of usages denoting high levels of engagement (time spent, number of actions carried out, diversity of types of documents consulted, etc.). Even though they are in a minority as far as overall numbers are concerned, they are far from being negligible in terms of absolute value, reflecting richly varied usages of the interface that would be useful to understand better. For example, although most Gallica sessions are single-type, 3% of sessions covering more than 5 documents go from one type to another: how are such clickstreams to be better understood? What actions and referring sites characterise them? What currently fosters such exploration of diversity?

Having a more refined semantic classification of documents available

For Gallica: apart from document types and other information contained in the OAI, and so far used separately (author, theme, etc., which are not always relevant or provided for all documents), it would be useful to be able to analyse research topics and their evolution over the course of clickstreams by having a record of the content of documents consulted available (e.g. via a number vector). An initial attempt at classification was made based on words in OAI notes; another method was also suggested ("word embedding", learned of from the corpus of articles in Wikipedia) but was not implemented due to lack of time.

For Data BnF: in order to gain a more detailed understanding of the makeup of Data BnF's audience and confirm or otherwise the "long tail" effect of consulting it, the fame of authors and/or works might be assessed in both conventional and automated fashion through use of information either internal (number of documents "about this author") or external (Wikipedia) to the site.

Incorporating the notion of user

Beyond the notion of "session", it would be useful to analyse user behaviour patterns with the help of a domain cookie with a short lifespan (1 month, for example). The same user might well make varied use of our interfaces, alternating short sessions and long sessions, etc. This possibility will be investigated in 2018.


Guillaume Ducourthial

Jean-Sébastien Affagard

Margaux Schmeltz

Xavier Solinas

Maeva Lopez-Poncelas

Christelle Bonod-Bidaud

Ruth Rubio-Amador

Florence Ruggiero

Jean-Marc Allain

Emmanuel Beaurepaire

Marie-Claire Schanne-Klein
email: marie-claire.schanne-klein@polytechnique.edu

Monitoring dynamic collagen reorganization

Keywords: multiphoton microscopy, polarization, biomechanics, collagen, skin G. Ducourthial et al.: Fast P-SHG imaging of collagen

Introduction

Collagen is the main component of the extracellular matrix and is organized as a network of fibrils in many tissues. [1,2] The distribution of this fibrillar collagen at the microscopic scale determines to a large extent the biophysical and biomechanical properties of tissues. However, while the mechanical behavior of tissues such as skin have been thoroughly characterized at macroscopic scale, the precise relationship between these macroscopic properties and the microstructure of collagen fibers is still a debate because of technical issues in acquiring multiscale experimental data. [3][4][5][6][7][8] In situ characterization of the collagen distribution in intact tissues is therefore a major issue to decipher the biomechanics of healthy and pathological tissues.

To address this issue, we have recently implemented a multiscale setup combining a traction device and Second Harmonic Generation (SHG) imaging. [9,10] SHG microscopy has indeed become the reference technique for 3D imaging of collagen at the microscopic scale. [11] SHG highlights collagen fibers with high specificity and sensitivity in intact tissues without requiring any staining. [12] Imaging of the collagen organization can be further improved by using polarization-resolved SHG microscopy (P-SHG). [13] This technique combines the sensitivity of polarimetry to molecular orientation and order with the specificity of SHG for fibrillar collagen. P-SHG measurements therefore provide additional information about the orientation and hierarchical structure of collagen fibers in the tissue. [14][15][16][17][18][19] However, P-SHG microscopy requires the acquisition of an image sequence recorded with a series of linear incident polarizations and is therefore a slow technique. Reliable P-SHG-based determination of the collagen reorganization during dynamic processes is therefore impossible if the sample exhibits distortions on the time scale needed to perform measurements with different polarizations. A number of strategies have been developed to overcome this limitation and have proven efficient in specific samples. [20][21][22][23][24][25] Nevertheless, most of these approaches are based on polarization analysis of forward SHG signals, which is not applicable to thick scattering tissues, or they use only two incident polarizations for SHG imaging, which limits their accuracy in tissues with disordered collagen. A robust method enabling fast epidetection of P-SHG signals is therefore needed to monitor collagen dynamic reorganization in thick deforming tissues.

In this paper, we describe an original fast P-SHG imaging setup based on line-to-line switching of polarization using an electro-optical modulator (EOM). As the EOM is located outside the SHG microscope, a careful calibration is performed to obtain reliable P-SHG measurements.

We use this P-SHG system to map the dynamic reorganization of the collagen network in ex vivo murine skin dermis during mechanical assays. We show that our fast setup enables monitoring collagen orientation even in this deforming sample and provides an accurate measurement of collagen reorganization during skin stretching.

Materials and methods

Fast P-SHG based on line-to-line EOM polarization switching

Multiphoton imaging is performed using a custom-built upright laser scanning microscope. [9] Excitation is provided by a femtosecond laser set at 860 nm (80MHz, 100fs, Mai-Tai, Spectra-Physics), with typically 20-30 mW excitation power at the focus of a water-immersion 20×, 0.95 NA objective lens (Olympus), resulting in 0.6 µm (lateral)  3.2 µm (axial) resolution measured near the sample surface. Multiphoton signals are epi-detected using photon-counting photomultiplier tubes (P25PC, Electron Tubes) and appropriate dichroic mirrors and spectral filters.

An Electro-Optic Modulator (EOM) (M350-210-02, Conoptics) is used to achieve fast rotation of the incident linear polarization as in the pioneering work by Stoller et al. [13] The incident polarization is set at 45° from the EOM axis using a half-waveplate (ACWP-700-1000-06-2, CVI) and the EOM is followed by a quarter-waveplate (ACWP-700-1000-06-4, CVI) oriented at 45° from the EOM axis (Figure 1.A). The resulting polarization orientation angle is then directly proportional to the input EOM voltage and can be tuned in the range 0°-170° at 860 nm. Polarization switching is achieved in less than 11 µs using a fast power supply (302A, Conoptics). It is synchronized with the microscope line scanning as displayed in Figure 1.B.

A quarter-waveplate is inserted at the back pupil of the objective lens to compensate for the ellipticity introduced by the microscope optical components. After optimization of all waveplate orientations, the ellipticity Emin/Emax is less than 10 % for all polarization angles and the power variation with the polarization angle is less than 1% (see Figure S1 in Supporting Information). Nevertheless, for any given setting, the scanning process causes a slight variation of the polarization angle across the field of view. This effect is calibrated by imaging a homogeneous fluorescent plate with a rotating polarizer inserted at the back pupil of the objective lens and the P-SHG images are corrected accordingly (see Supplementary method and Figure S2 in Supporting Information).

Alternatively, slow P-SHG imaging can be performed by use of a motorized half-waveplate at the back pupil of the objective as reported previously. [26] In this case, as mechanical rotation is slow, polarization orientation changes are performed between two successive images instead of between lines.

Traction assays in murine skin under P-SHG microscope

All animal experiments are performed under animal care procedures in accordance with the guidelines set by the European Community Council Directives (86/609/EEC). All experimental procedures are approved by the Direction of the Veterinary Service of Rhone Department (DDSV, Lyon, France).

Reorganization of the collagen network is monitored in ex vivo murine skin dermis during traction assays as previously described. [9,27] Eight wild-type mice (129sv) are sacrificed by cervical dislocation at one month. Skin from the back is depilated, the epidermis is removed using 3.8 % ammonium thiocyanate [28] and the remaining dermis is cut into a dog-bone shape along the antero-posterior direction to ensure homogeneous uniaxial tensile load in the central tested portion. [9] Skin samples are attached to a custom-built uniaxial traction device using metallic jaws and inserted in place of the microscope stage, with the papillary dermis up facing the objective lens. Symmetric stretching is applied at a slow strain rate (typically traction speed of 2 µm.s -1 ) to monitor the same region of interest (ROI) during all the experiment while measuring the force every second. Stretch is maintained constant during multiphoton imaging, resulting in an incremental loading path with steps of 0.05 stretch ratio. No preconditioning is applied before stretching. Immersion gel (Lacrygel, Europhta) is used to ensure optical contact with the objective lens and to prevent skin dehydration during experiments. 375 × 375 × 50 µm 3 image stacks are epi-collected at 100 kHz pixel rate, using 0.5 µm pixel size and 2 µm axial steps. 18 different incident polarizations (0-170° every 10°) are used for P-SHG. The total acquisition time of all polarizations for a line of 750 pixels using EOM switching is 200 ms. Mechanical data are processed as previously described. [9,27] Briefly, the global stretch ratio λ is obtained as the imposed length between the jaws divided by the initial length, and the nominal stress as the measured force divided by the initial skin section. The tangent modulus is then obtained as the slope of the linear part of the nominal stress/stretch curve (see Supporting Information).

P-SHG theory and image processing

The SHG response of collagen fibers is described by a second order nonlinear susceptibility tensor . Considering that the collagen fibers exhibit a cylindrical symmetry and that the Kleinman symmetry applies, [13] there are only 2 independent components: and where x is the fiber axis. Derivation of the P-SHG signal intensity from a collagen fiber oriented at an angle φ in the image plane XY then requires to change from the fiber frame xyz to the microscope reference frame XYZ and to calculate the SHG response to a linearly polarized excitation at an angle θ to the X axis (see inset in Figure 1(A)). It writes in the plane wave approximation, without any polarization analysis of the SHG signal because of the epidetection configuration:

| |² sin 2 ² (1)
Equation ( 1) also applies to fibers lying out of the image plane and projected in the image plane at an angle φ to the X axis, since the out of plane dependence is included in the susceptibility components in the XYZ frame. The parameter K includes various geometrical parameters and the square of the excitation laser intensity. Note that this analytical approach does not take into account the axial components of the excitation field due to strong focusing. [29] Indeed, we have recently demonstrated that strong focusing does not affect the determination of the in-plane orientation φ of collagen fibers, while it strongly affect other parameters, notably the ratio of the two susceptibility components. [26] It is convenient to write Eq. ( 1) as a function of its Fourier components: [13,30]

θ 2 4 (2)
Where A, B and C are related to the susceptibility components of the collagen fiber. In this study, we are only interested in the in-plane orientation φ of the collagen fiber. We use a FFT algorithm to determine the A, B and C complex parameters in every pixel and the orientation φ is obtained from the phases of B and C. [26] The reliability of this determination is assessed by a coefficient of determination, so-called R 2 parameter, that compares the experimental data and the curve obtained from the A, B and C parameters extracted with FFT. Angles determined with R 2 <0.5 are skipped in the quantitative analyses, which ensures an angular accuracy of 5°. This accuracy is estimated by the difference between the first angle φ extracted from B and the second angle φ extracted from C. The orientation maps are displayed using a HSV Look-Up-Table : the Hue (H) indicates the collagen orientation φ and the brightness (V) the R 2 parameter:

V=1 for R 2 ≥ 0.5 and V [0;1] for R 2 [0;0.5].

Analysis of the reorganization of the collagen fiber network

The resulting map of fiber orientations is further processed to plot a histogram of the orientations within the field of view. Our previous studies have shown that tissue relaxation is negligible during the acquisition of the P-SHG image, [27] so that we can safely quantify parameters from this histogram, namely calculate the entropy S of this orientation distribution.

S is defined as the usual statistical entropy, using angular bins of 1° or 10°, and is normalized to be independent of the number of bins nb(θ):

∑ ln ° ° (3)
Where p(θ) is the normalized number of pixels with the orientation θ in the image, which is directly obtained from the orientation histogram. The entropy S then assesses the degree of orientation disorder in the field of view: S=1 for a fully isotropic orientation distribution and decreases towards 0 for thinner orientation distribution. This analysis is performed for the six consecutive 2D images exhibiting the maximal SHG intensity in the P-SHG z-stack. The mean ± SD (standard deviation) entropy is calculated and plotted as a function of the stretch ratio for every skin sample. Linear fitting gives the slope of the entropy variation that characterizes the reorganization of the collagen fiber network upon stretching.

Alternatively, the entropy is computed from orientation maps obtained by image processing of the SHG images using the morphological filtering method reported previously. [9] Briefly, all the SHG images obtained for different linear excitations are summed up to mitigate the impact of polarization versus fiber orientation on the SHG intensity. The local orientation of the fibers is then extracted in every pixel using morphological filtering by a rotating linear structuring element (19 pixels/9.5 µm long). This method provides histograms with 10° bins. [9] 3. Results

Fast P-SHG imaging allows orientation measurements in dynamic samples

P-SHG imaging was performed in the dermis of eight wild-type mice during traction assays and the collagen orientation was processed pixel-wise using Eq. ( 2). Note that faster motion artefacts in the observed tissue would not be removed by this method and would necessitate decreasing the number of pixels in order to decrease the acquisition time.

Nevertheless, such fast and strong movements are only observed sparsely at the very beginning of the traction assay, probably due to unfolding of skin wrinkles, which means that they do not affect the linear part of the collagen reorganization (see next section).

Fast P-SHG imaging of collagen fiber network reorganization

Fast P-SHG was used to assess the dynamic reorganization of collagen fibers in ex vivo murine dermis using sequential 0.05 stretching steps. Force measurements were processed to plot the stress-stretch curve and derive the tangent modulus as shown in Supporting Information (Figure S3). We obtained a J-shaped curve as in our previous studies, [9] which means that the increased duration of the traction assays due to polarization resolution did not strongly affect the mechanical response of the dermis at macroscopic scale. The tangent modulus is slightly larger in this study: 2.1 ± 0.2 MPa (N=8) than in the previous one: 1.6 ± 0.1 MPa (N=31), but the difference cannot be considered as statistically significant (p=0.07 using Wilcoxon test). [10] Figure 3.A displays typical sum-SHG images at a few stretch ratios, the loading axis being the horizontal direction. All these images display a similar pattern of elliptical features without any SHG signal, which correspond to the positions of the hair follicles. It shows that we successfully imaged the same region of interest at all stretch ratios. The deformation of this pattern upon stretching enables the computation of a local stretch ratio, in excellent agreement with the imposed global one (data not shown). [9] It shows that the dermis response is homogeneous and validates the mechanical assays.

Figure 3.A also shows that the collagen fibers gradually align along the traction direction at increasing stretch ratios. This effect is better revealed in Figure 3.B, which displays the orientation maps obtained from the P-SHG images. The maps at low stretch ratios display various colors corresponding to various orientations, while the maps at high stretch ratios gradually turn red that is the color of the horizontal traction direction. which corresponds to the antero-posterior direction, even at low stretch ratio, as already reported. [27] The observed peak is quite wide at low stretch ratios and becomes gradually thinner at higher stretch ratios. Accordingly, the entropy of these orientation distributions decreases when the stretch ratio increases, as shown in Figure 4.B. Note that the entropy is sufficient to fully characterize the distribution here as the mean orientation of the collagen fibers is aligned with the stretching direction. We observe a linear decrease, which is preceded by a plateau at low stretch ratios in some skin samples. This observation is in agreement with our previous experiments that used morphological filtering to extract orientation maps from SHG images acquired with circular polarization. [9] It indicates that the increased duration of the experiment due to polarization resolution does not strongly affect the reorganization of the collagen microstructure in murine skin samples. The slope obtained by linear fitting then provides a quantitative measure of this reorganization, as summarized in Table 1 for the eight skin samples.

Comparison with morphological filtering

For the sake of comparison, we also processed the sum-SHG images using morphological filtering as previously reported. [9] The corresponding orientation maps, orientation histograms and entropy evolution as a function of the stretch ratio are depicted in the Figure S4 (Supporting Information), and the resulting slopes for the eight murine samples are summarized in Table 1.

These orientation results show the same trends as the ones obtained from P-SHG, but the orientation maps derived from the sum-SHG images show more different orientations and the histograms are wider and exhibit higher background. As a consequence, the values of the entropy and of the slope of its evolution upon stretching are different from the ones obtained from P-SHG. [START_REF] Serra | Analysis and Mathematical Morphology[END_REF][START_REF] Quinn | [END_REF] We also compared the respective precision of P-SHG and morphological filtering by looking at the accuracy of the linear fitting of the entropy decrease upon stretching. This accuracy is assessed by calculating the R² parameter of the linear fit. It is slightly larger for P-SHG (0.95 ± 0.02) than for morphological filtering (0.88 ± 0.03).

Discussion and conclusion

In this study, we implement a fast P-SHG microscope based on an EOM, following the pioneering work by Stoller et al. [13] We take advantage of the fast switching rate of the EOM to synchronize the rotation of incident linear polarization with line scanning and record XPY images. Our microscope advantageously enables imaging of thick scattering samples because it is based on epidetection without polarization analysis. It may easily be generalized to circular polarizations as well. Note however that a careful calibration of the polarization states is necessary because the EOM device is conveniently installed outside of the microscope, resulting in the propagation of the incident polarization states across potentially distorting components of the microscope. Nevertheless, once the calibration is completed, this fast P-SHG microscope is a robust and application-friendly setup.

Importantly, the acquisition duration depends on several parameters here: the number of pixels per line, the pixel dwell time and the number of incident polarization states. While the total acquisition time of a XPY image is longer than with circular-polarization (150 s versus 8.4 s for 750 x 750 pixels), the relevant time here is the one for line imaging: 200 ms for 750 pixels.

This acquisition duration can be easily minimized by decreasing the number of pixels in every line, either using larger pixels (faster scanning) or imaging smaller areas. For instance, lines of 100 pixels instead of 750 pixels would result in a P-SHG acquisition time of 27 ms per line instead of 200 ms. This acquisition time is only moderately longer than the one relevant for morphological filtering in one pixel, that is the acquisition time of the window used for filtering (in our case, 19 x 19 pixels): 5.4 ms. Some studies have proposed to further reduce the acquisition time by using only 2 or 3 polarization states either for the excitation or for the analysis of circularly-excited SHG. [21,22,24] While this approach may be efficient in samples with simple orientation distribution, it cannot provide accurate orientation maps in samples with more complex orientation distribution.

In practice, the pixel dwell time and the number of incident polarization states cannot be easily decreased because they determine the signal to noise ratio in the P-SHG images. Indeed, the accuracy of the determination of the collagen orientation and of other parameters extracted from P-SHG images is related to the total number N of detected photons. [33,34] This number N scales as the pixel dwell time multiplied by the number of incident polarization states. It varies in every biological tissue: it mainly depends on the maximum intensity usable for imaging without photo-damaging the tissue and on the collagen density and alignment in the tissue, which determines the efficiency of the SHG process. It is also smaller in depth because scattering and aberrations may decrease the excitation intensity. Practically, in our murine skin dermis samples, we safely use 20-to-30 mW excitation power and detect approximately N=500 photons in total per pixel. This corresponds to an accuracy of 3° on the orientation determination, which is slightly degraded to 5° (see R² threshold in section 2.2) by the residual ellipticity of the excitation and mostly by possible polarization distortions due to heterogeneities within the skin sample. [26] Compromises on this accuracy would enable to speed up the acquisition, for instance 7 ms for one line of 100 pixels using 9 polarizations instead of 18 ones and half pixel dwell time (5 µs). Nevertheless, given the order of magnitude of SHG signals detected in collagen tissues, the minimal pixel dwell time in a point-scanning microscope is a few µs. It means that in collagenous tissues, the main limitation to decrease the acquisition duration is the relatively low SHG signal, not the polarization switching rate. Smart approaches to achieve video-rate polarization switching have been reported recently and proved efficient in samples exhibiting higher signals. [23,25,35] Nevertheless, in complex biological samples such as skin dermis, such imaging rates require to rely on subsequent frame averaging to achieve a sufficient signal to noise ratio and result in a similar acquisition duration as the one we obtained for full images. In that respect, our strategy based on fast line-to-line P-SHG imaging is a good compromise: it is slow enough to detect small endogenous SHG signals and fast enough to determine collagen orientation in dynamic samples, as demonstrated in Figure 2. Moreover, it works in thick samples where epi-detection is mandatory, which is not the case for many other approaches. [22][23][24][25]35] It may be further improved by synchronizing the EOM with pixel imaging in order to image tissues with faster motion. In case of even faster tissue deformations, processing SHG images acquired with circular polarization, short pixel dwell time and low sampling may be the only possible technique to determine collagen orientation, probably at the expense of the accuracy.

We use here this fast P-SHG microscope to characterize the reorganization of the collagen network in ex vivo murine dermis upon stretching. We previously performed similar experiments using morphological filtering to determine the collagen orientation from SHG images, and our results in wild-type and genetically-modified mice challenged the usual interpretation of the microstructural origin of the skin macroscopic biomechanical behavior. [27] Indeed, as we have previously reported [9,10,27], the collagen fibers do not align with the direction of traction in the non-linear part, but they do so in the later linear part of the stress/stretch curve. These observations are in contradiction with the classical assumption of the role of collagen in skin mechanical properties. [3][4][5] Nevertheless, advanced modeling of the collagen network response requires measurements of the collagen reorganization with a good spatial accuracy that is pixel-wise measurements, while methods based on image processing probe a larger scale. Indeed, our method based on morphological filtering of sum-SHG images provides a pixel-wise orientation map, but it probes in practice the collagen orientation at the scale of the rotating linear structuring element (here 19 pixels, that is approximately 10 µm). It is the case for all image processing techniques aiming at determining orientations, which always require to process a window around the pixel under study. [START_REF] Serra | Analysis and Mathematical Morphology[END_REF][START_REF] Quinn | [END_REF] In contrast, P-SHG effectively probes the collagen orientation at the pixel size (here 0.5 µm).

P-SHG measurements are therefore required to improve the spatial accuracy of collagen orientation measurements during traction assays in the framework of our multiscale study.

Accordingly, we successfully image the reorganization of the collagen network in murine dermis using fast P-SHG imaging. We obtain thin orientation histograms with low background that show qualitative similarities to the wide histograms with high background obtained by morphological filtering. The observed quantitative differences are attributed to the different scales probed by the two methods. Moreover, both methods have to compromise between (i) the accuracy of the orientation maps (number or angular bins, pixel size and binning, R² filtering…) and (ii) the number of pixels with determined orientation and consequently the statistical significance of these orientation maps. This compromise does not imply the same parameters in both methods and is therefore difficult to adjust in the same manner, while the set of parameters used in every method may change the quantitation of the results. For instance, the results obtained from P-SHG images using 1° and 10 ° angular bins for the calculation of the entropy are quantitatively different. Nevertheless, the results obtained by P-SHG (using 1° or 10° bins) and by morphological filtering reveal the same trends, as shown by the last columns of Table 1: the ratio of the entropy slopes obtained by both methods exhibits low variation for the eight skin samples. Notably, these entropy slopes obtained by P-SHG are determined with a better accuracy than the ones obtained by morphological filtering. Note however that this conclusion holds for dermis, a tissue with a dense network of collagen fibers, but the situation may be different in other types of tissues, for instance in the case of sparsely distributed straight collagen fibers, which orientation is accurately determined by SHG image processing. [START_REF] Quinn | [END_REF]36] These results are therefore a proof of feasibility of fast P-SHG imaging of skin microstructure reorganization during mechanical assays. Systematic analysis of more skin samples should therefore enable the determination of accurate orientation histograms and enable detailed comparison with multiscale mechanical models.

In conclusion, we report here on a new fast P-SHG microscope based on EOM polarization switching synchronized with line imaging. This method is fast enough to obtain reliable orientations in dynamic collagen samples, while slow enough to get a good SHG signal to noise ratio in thick collagenous tissues. We demonstrate that it enables accurate quantitation of dynamic collagen reorganization in murine skin dermis during stretching biomechanical assays.

This approach can in principle be generalized to other multiphoton modes of contrast and enable fast imaging of a variety of orientation-sensitive processes. (same as in Figure 3). A: Orientation histograms obtained from P-SHG orientation maps using 1° bins at increasing stretch ratios. The baseline is shifted upwards for increasing stretch ratios.

B: Variation of the normalized entropy of the orientation distributions in A as a function of the stretch ratio. The black squares correspond to the experimental data: mean ± SD calculated on the 6 2D P-SHG images exhibiting the maximal signal; the solid line corresponds to linear fitting.

Table 1. Measurements of dermal collagen reorganization upon stretching for 8 wild-type mice using morphological analysis and P-SHG method. The slope of the entropy decreases upon stretching, and the R² of the linear fit are provided for every mouse using both methods, either for 10° angular bins (morphological filtering and P-SHG) or for 1° angular bins (only P-SHG).

The mean  SEM of every parameter is also calculated.

Graphical Abstract

An original fast polarization-resolved SHG microscope based on line-to-line switching of polarization using an electro-optical modulator is implemented to map the dynamic reorganization of the collagen network during mechanical assays. This system successfully highlights the collagen alignment along the traction direction in ex vivo murine skin dermis and enables quantitation of the linear increase of the collagen orientation distribution entropy as a function of the stretch ratio.

Figure 2

 2 displays the orientation map for the WT4 mouse at 1.05 stretch ratio, as well as the SHG image obtained as the summation of all images acquired with the set of different linear incident polarizations (sum-SHG image, Figure2.A). This sum-SHG image allows the visualization of the collagen fibers regardless of their orientation in the image plane, in a similar way as the SHG image obtained by use of circularly-polarized excitation (circular-SHG image, Figure2.C). Nevertheless, the circular-SHG image is recorded much faster (8.4 s for 750x750 pixels at 100 kHz pixel rate) than the sum-SHG image that requires the acquisition of 18 polarizations for every line (18 x 8.4 s = 150 s). Accordingly, motion artefacts are observed in the sum-SHG image, not in the circular-SHG one. Nevertheless, these line-to-line motion artefacts do not impede the determination of the collagen orientation because P-SHG images are acquired using fast rotation of polarization synchronized with line scanning. The same line is recorded sequentially with 18 different incident linear polarizations using fast polarization switching by the EOM device, then the scanner moves to the next line. This method enables fast recording of P-SHG lines (200 ms), even if the full image is recorded much slower (200ms x 750 lines = 150s, see above). Accordingly, Figure2.B shows that the collagen orientation is satisfactorily processed from the P-SHG images even in the regions with motion artefacts.

Figure 4 .

 4 Figure 4.A displays the orientation histograms obtained from these orientation maps at increasing stretch ratios. The collagen fibers show a main orientation close to the traction axis,

Figure 1 .

 1 Figure 1. Experimental setup. (A) Fast P-SHG microscope using an EOM and a quarterwaveplate to rotate the incident linear polarization. The imaging plane geometry is shown in the insert. (B) Timing diagram. EOM voltage is synchronized with line scanning (x direction). The next line is scanned (y scanning) when the first line has been scanned for all EOM voltages (i.e. all polarizations).

Figure 2 .

 2 Figure 2. Zoomed-in images of a region with sample motion (ex vivo skin dermis of WT4 mouse at 1.05 stretch ratio). A: summation of all P-SHG images acquired at various linear incident polarizations; B: processed orientation image obtained from P-SHG acquisition; C: SHG image acquired by use of circular polarization. Note that the last image is slightly different from the first one because it was acquired afterwards and collagen reorganized. Scale bar: 20 µm.

Figure 3 .

 3 Figure 3. P-SHG imaging at increasing stretch ratio: 1.1, 1.2, 1.3 and 1.4 in skin sample WT15. A: Sum-SHG images obtained as the summation of all P-SHG images obtained for the different linear incident polarizations. B: Orientation maps processed from the P-SHG images and plotted using the HSV look-up table: H corresponds to the orientation shown in the inset, and V corresponds to R 2 : V=1 for R 2 ≥ 0.5 and V[0;1] for R 2 [0;0.5]. Scale bars: 50 µm.

Figure 4 .

 4 Figure 4. Reorganization of the collagen fibers network upon stretching in skin sample WT15

Acknowledgements

The authors acknowledge financial support from Agence Nationale de la Recherche under contracts ANR-10-INBS-04 (France BioImaging), ANR-11-EQPX-0029 (Morphoscope2) and ANR-13-BS09-0004-02 (Metis).

Supporting Information

Additional supporting information may be found in the online version of this article at the publisher's website. 3 and4). The nominal stress (measured force divided by the cross-section) is plotted as a function of the stretch ratio (imposed length divided by the initial length of the sample). The tangent modulus is obtained as the slope of the linear part.

Figure S3. Mechanical response of skin sample WT15 (same as in Figures

Figure S4. Reorganization of the collagen fibers network upon stretching in skin sample

WT15 (same as in Figures 3 and4). A: Orientation maps at 1.1, 1.2, 1.3 and 1.4 stretch ratios processed from the sum-SHG images using morphological filtering. Look-up table: HSV, where H codes for the orientation and V for the intensity below Imax/2 (V=1 for I>Imax/2, V[0;1] for I[0;Imax/2]). Scale bar: 50 µm. B. Orientation histograms obtained from morphological filtering using 10° bins at increasing stretch ratios. The baseline is shifted upwards for increasing stretch ratios. C: Variation of the normalized entropy of the orientation distributions in B as a function of the stretch ratio. The black squares corresponds to the experimental data: mean ± SD calculated on the six 2D SHG images exhibiting the maximal signal; the solid line corresponds to linear fitting.


Simon Gorecki
email: simon.gorecki@u-bordeaux.fr

Judicaël Ribault
email: jribault@gmail.com

Grégory Zacharewicz

Nicolas Perry

Yves Ducq
email: yves.ducq@u-bordeaux.fr

Judicael Ribault

Gregory Zacharewicz
email: gregory.zacharewicz@mines-ales.frperry.nicolas@ensam.eu

Distributed simulation for a modeling and simulation tool: PAPYRUS

Keywords: Model Driven Architecture, HLA, FMI, Papyrus, Simulation

come

INTRODUCTION

In industry, the design of complex systems requires Modeling & Simulation (M&S) steps to represent behaviors and interactions in the system and between systems. As technologies are growing, the systems complexity also increases and makes systems more difficult to model and simulate. Along with this growing complexity, risks, hazards, and threat must also be considered during the modeling process. Several researches which prove this fact have been done about project risk management [START_REF] Altuhhova | An Extension of BPMN for Security Risk Management[END_REF][START_REF] Better | Simulation Optimization Applications in Risk Management[END_REF]. In this context, it is required to divide the complexity of a system using Distributed Simulations (DS). However, industrial context providing the application case is using UML and Papyrus that are not ready for DS.

Papyrus, an open source UML / SysML modeler, provides to users and developers a tool to specify UML models and then to execute them thanks to the Foundation UML Subset (fUML) standard (OMG Publication 2018). Papyrus is already used by our partner company to model and simulate the production of solar power plant. This software is used thanks to its user friendly modeling and simulation features according to the fUML standard. We can also describe specific behaviors according to the UML profile. Those profiles can be taken into account during the simulation run thanks to the Moka execution engine. However, as introduced before, growing complexity implies many risk and hazard impacts on models. Modeling all these risks would be very complex and make the initial model unreadable [START_REF] Gorecki | Including Co-Simulation in Modeling and Simulation Tool for Supporting Risk Management in Industrial Context[END_REF].

In order to solve this issue, [START_REF] Gorecki | Including Co-Simulation in Modeling and Simulation Tool for Supporting Risk Management in Industrial Context[END_REF] first built an extension of the Moka execution engine to handle risk and hazard outside the model with a Functional Mock-up Unit (FMU). The FMI standard allowed us to relocate and factor complex tasks without overloading models. In the following of these works, to deal with risks management in M&S, we proposed to divide complexity into several DS components. The purpose of this paper is managing Papyrus instances as DS components to allow users and developers to deal with complexity and time constraints. In details, after having managed interaction with the Moka scheduler in previous works, the goal is now to implement synchronization points mechanism in Moka. This is a complex problem since it implies defining a synchronization, communication algorithm which is not defined by the FMI Co-Simulation standard.

To solve the time related issues in distributed context, we use the High Level Architecture (HLA) standard (IEEE Computer Society 2010a, 2010b, 2010b), one of the most used standard in DS. Indeed, in HLA, the Run Time Infrastructure (RTI) can solve interoperability issues between components, provide communication mechanism and time management. Hence, to create a Papyrus distributed simulation, we propose to link HLA and the Moka engine through FMI. We have already coupled Moka and FMI in our previous paper, and some research already exist on coupling FMI and HLA [START_REF] Bouanan | Including in HLA Federation Functional Mockup Units for Supporting Interoperability and Reusability in Distributed Simulation[END_REF][START_REF] Sievert | Adapting Functional Mockup Units for HLA-Compliant Distributed Simulation[END_REF]. In short, the main purpose of this contribution consist in merging these two standards, HLA RTI as scheduler for several Papyrus instances and FMI for compatibility with other components.

BACKGROUND

This section introduces the different concepts used in this paper for the proposed solution. The first subsection is dedicated to the Papyrus tool and its execution engine: Moka. Then, the Functional Mock-up Interface standard (FMI) is described. Finally the High Level Architecture (HLA) standard is explained in 2.3 section.

Papyrus & Moka execution engine

Papyrus is the UML/SysML modeler of the Eclipse foundation. UML allows the user several modeling notations to represent structures and behaviors of its system. This standard is however not very efficient for system description. That is the reason why the UML profile exists. Created by the Object Management Group (OMG), a UML profile allows users to specify its own UML language on its own semantic objects.

Papyrus provides tools for modeling, and executing these two standards. It is based on Eclipse and is an open source project. The execution part is handled by MOKA, a fUML execution engine. In accordance with its primary goal to implement the complete standard specification of the UML2, Papyrus provides an extensive support for UML profiles.

Moka is a Papyrus plugin designed to provide an execution environment complying with the OMG standards UML. The Foundational UML Subset (fUML) is an executable subset of the UML standard which can be used to define, in an operational style, the structural and behavioral semantics of systems in order to be simulated. [START_REF] Guermazi | Executable Modeling with FUML and Alf in Papyrus: Tooling and Experiments[END_REF]. With this standard, Moka is able to execute UML models designed under the Papyrus UML editor. Moreover, designed for Papyrus, Moka is open source and can be extended to support the UML Profile or alternative execution semantics, and thereby, be adapted to multiple usages.

Functional Mockup Interface (FMI)

The Functional Mock-up Interface (FMI) for Co-Simulation interface is designed to simulate tools coupling and subsystem models coupling. In both cases, coupled systems are exported by their own tools. FMI is a standard interface for time dependent coupled systems which are continuous in time or time-discrete [START_REF] Bastian | Master for Co-Simulation Using FMI[END_REF][START_REF] Blochwitz | Functional Mock-up Interface for Model Exchange and Co-Simulation[END_REF][START_REF] Sievert | Adapting Functional Mockup Units for HLA-Compliant Distributed Simulation[END_REF]. It provides interfaces between master and slaves and addresses both data exchange and algorithmic issues. FMI for Co-Simulation consists of two parts:

 Co-Simulation Interface: a set of C functions for controlling the slaves and for data exchange of input and output values as well as status information.  Co-Simulation Description Schema: defines the structure and content of an XML file. This slave specific XML file contains "static" information about the model (input and output variables, parameters …) and the solver/simulator (capabilities …). The capability flags in the XML file characterize the ability of the slave to support advanced master algorithms which use variable communication step sizes, higher order signal extrapolation etc.

A component implementing the FMI is called Functional Mock-up Unit (FMU). It consists in one zip file containing the XML description file and the implementation in source or binary form (dynamic library). A master can import an FMU by first reading the model description XML file contained in the zip file. Coupling simulators by FMI for Co-Simulation hides their implementation details and thus can protect intellectual property.

High Level Architecture

HLA defines a framework which allows the creation of a Distributed Simulation (DS). DS participants are called federates, they can communicate with each other. The standard was originally created by the Office of Defense Modeling and Simulation (DMSO) of US Department of Defense (DoD) to facilitate the assembly of stand-alone simulations with a different architecture. The original goal was the reuse of military applications, simulations, and sensors. This standard is designed to resolve interoperability issues and reusability between software components. Another interesting aspect of this specification is the synchronization aspect. It allows the dynamic management of interoperability issues with simulations exchange messages: it must be ensured that messages are sent at the right time, in the right order, and that they do not violate causal constraints. To do so, various systems for synchronization of processes and time management are proposed by HLA.

According to the HLA standard, all simulations participating to the application are called "federate". A classical HLA federate consists of a simulation model and local RTI component (LRC). The simulation model is a physical, mathematical, or logical representation of processes and systems. These entities can communicate with each other through a Run-Time Infrastructure (RTI). It is the RTI which manages the federation, authorizes federates to communicate or not, and provides various services such as time management, file or data exchange, etc. The FOM file is a XML file which describes interactions/communications between federates.

STATE OF ART

HLA and FMI are well known standards of Co-Simulation and Distributed Simulation. They both are described in different specification books. For HLA, (IEEE Computer Society 2010a) describes how interactions can be done between each component of de DS: Exchange mechanism are done by Subscribe (read) and Publish (write) provided by the RTI. Those communications are described by SOM and FOM XML files. However, there is no simple way to describe the federate behavior inside the federation: each federate is autonomous.

HLA imposes an architectural rule for its federates. Each one of them has to assume its own time management. HLA architecture can handle three advancing time mechanisms for each federate member of the federation. Authors presented in [START_REF] Bouanan | Including in HLA Federation Functional Mockup Units for Supporting Interoperability and Reusability in Distributed Simulation[END_REF] an application running several federates connected to the same RTI and communicating together. Among these federates, one of them is a hybrid federate hosting an FMU making a communication bridge between HLA classical federates and a hybrid HLA FMI element.

In literature, efforts has been done to upgrade the FMI standard by adding extensions and enable treatment of discrete-time models [START_REF] Franke | Discrete-Time Models for Control Applications with FMI[END_REF]. In this paper, authors proposed a solution to synchronize in time FMUs with others elements by adding a clock in the model description and implement interface to call it. An FMI export with synchronous features was implemented into the tools Dymola and OpenModelica.

An FMU import was implemented in the optimization solver HQP.

Other papers are also interesting concerning use of Papyrus software. [START_REF] Guermazi | Executable Modeling with FUML and Alf in Papyrus: Tooling and Experiments[END_REF] discuss about the fUML and the Alf standards available in Papyrus which allow to execute UML models. Authors of this paper highlight four major concerns regarding the design of fUML execution engine: Extensibility, control and observability, time support, and connectivity with external tools.

In [START_REF] Ellner | A FUML-Based Distributed Execution Machine for Enacting Software Process Models[END_REF], authors highlight the limits of the fUML standard by explaining is lack of support for distributed execution in order to guide and support team members with their activities. FUML also requires explicit modeling of auxiliary user specific attributes and behavior of model elements, which is a repetitive and error-prone task. The author presented in this paper fUML extensions able to support distributed execution and human interactions. With these FUML extensions it becomes feasible to enact reusable standard software process models in realistic projects.

In [START_REF] Gorecki | Including Co-Simulation in Modeling and Simulation Tool for Supporting Risk Management in Industrial Context[END_REF], authors proposed an entry point in the Moka execution engine allowing control on local scheduler. This paper demonstrates how to interact with time management during an execution process in order to simulate hazards. This Moka extension allows to outsource risk management in an independent tool in order to connect an FMU during the simulation execution process and outsource tasks. This concept will be more explained in the next section.

CONTRIBUTION

The contribution in this paper is based on two previous works recalled in state of the art. The first one is about controlling simulation time in Papyrus with a specific Moka extension interacting with an FMU file [START_REF] Gorecki | Including Co-Simulation in Modeling and Simulation Tool for Supporting Risk Management in Industrial Context[END_REF]). The second one is about creating an HLA federate hosting a Hybrid HLA/FMI federate able to communicate with each other [START_REF] Bouanan | Including in HLA Federation Functional Mockup Units for Supporting Interoperability and Reusability in Distributed Simulation[END_REF]. Merging these two contributions consists in defining a new Moka extension able to communicate with another instance of Papyrus, allowing communication and time synchronization through FMI and HLA mechanism. Section 4.1 is recalling previous works to situate the contribution of this paper, while the contribution is detailed in section 4.2 and 4.3.

Related contributions

The main objective of this paper is to make a Papyrus simulation able to communicate with others Papyrus instances at runtime. This contribution is based on a previous paper [START_REF] Gorecki | Including Co-Simulation in Modeling and Simulation Tool for Supporting Risk Management in Industrial Context[END_REF] where we proposed a Moka extension able to communicate with an external tool (see Figure 1). This tool had to generate hazards and insert them into the Papyrus simulation process. This previous paper demonstrated that we can use an external tool to influence Papyrus simulation time. The Moka extension has access to the content of our current custom UML profile, as well as the whole simulation. In this extension, we are able to define (with java code) specific behaviors according to what is given by the UML profile. In Figure 1 above, the Moka extension makes a request to an external FMU. This FMU can generate hazards and delay tasks depending on parameters which are the content of the UML profile. Figure 2 describes a simple example of how the UML profile and the Moka extension can affect a simple model. The left sequence flow diagram shows the normal execution of Figure 1 without applying our Moka extension: each task is executed according to its "normal" time (normal time is the time affected to a task in fUML Papyrus model). In the right side of the figure, we applied our UML profile on Task 3 (as in Figure 1). We can see when Task 3 is executed, the Moka extension sends all the information contained in the UML profile to the FMU. The FMU determines (based on the information contained in the UML profile) if this task should be delayed or not, and if yes, for how long. We can see on the right sequence flow diagram that the Moka extension extends the simulation time of Task 3 to 3 units of time. In this example, the FMU determine delay time according to content of the UML profile (name of the task, simulation time, etc) and use databases, web services, and statistics algorithme to determine the delay. [START_REF] Bouanan | Including in HLA Federation Functional Mockup Units for Supporting Interoperability and Reusability in Distributed Simulation[END_REF]) presents a hybrid HLA/FMU connected to an RTI, which communicates with an external application through the FMI standard. In our paper, we use the same approach: a hybrid HLA/FMU component able to communicate with one local Papyrus execution and a HLA federation.

Hybrid HLA/FMI federates to manage Papyrus instances

The objective here is to run several Papyrus instances, each one having its own hybrid HLA/FMU component. We use the FMI standard as entry point in the Papyrus simulation, and the HLA standard for its functionality of time management, and communication protocols. The main objective is to implement waiting points between Papyrus instances. Using a hybrid component brings us several advantages. First, all the development which we are doing are be compressed in one single FMU file. Indeed, it is easy to reuse our proposition in futur Papyrus project directly on importing the file. Also, by assuming the FMI standard, we can connect other compliant tools in industry. And finally, using the HLA standard enable connection of other simulations to the final DS.

When an instance of Papyrus has to wait for another instance of Papyrus, messages are exchanged amongst them through HLA communication mechanism. At the beginning of the federation, each federate subscribes to an interaction class "Message" (described in Figure 2) in order to be informed of every state modification during the global distributed simulation. HLA mechanism does not allow a federate to address message directly to a recipient. Federates can only publish and subscribe Objects or Interactions. At the initialization of the federation, each federate subscribes to the same interaction class (Figure 4). Indeed, during the simulation, each federate receives all interactions. They use "src" and "dst" attributes to identify if there are concerned by the message or not. If not, they simply ignore the message. If yes, they can publish an answer. When an instance of Papyrus enter in a synchronized state, the associated federate publishes an interaction in order to inform the other members.

Each Papyrus instance has the same UML Profile and the Moka extension corresponding to the Figure 3 below.  idTask: a unique identifier. It will allow any federate to identify where it comes from. For example, P1.3 is the id of the Task 3 in the federate 1.  idLink: id of the associated locker/unlocker task.  isLocked: a Boolean variable that determine if it is a locking task or an unlocking task.

At the execution of a synchronized task, the Moka extension shares with the hybrid HLA/FMI federate information contained in the profile. The federate behaviors depend on the status of the system.  src contains the unique id of the message transmitter. If it is a "lock" or "unlock" message, it contains something such as "P1.3". If it is a reply to a lock message, it only contains its Papyrus id: "P1".  dst contains the target id of the message. With the HLA subscription and publication mechanism, it is impossible to address a message to a specific person. So the federate is able to filter which messages are for them and which are not (if not, ignore them).  localTime is given by the Moka extension. When a lock message is received, the recipient must reply by its own local time.  isLocked is about the semantic of the message. 0 corresponds to a unlocking message. 1 corresponds to a locking message. 2 corresponds to a locking reply (the federate has to give its Papyrus local time).

During the simulation, each federate has the opportunity to send:

 a "lock message" when Papyrus enters into a synchronized point which needs an unlocking from another Papyrus process.  an "unlock message" when Papyrus enters into a synchronized point which unlocks another Papyrus process.  a "time reply message" when a federate receives a notification from a locked Papyrus process, it has to reply its actual local time.

From a technical point of view, we are using in this application the HLA Run Time Infrastructure (RTI) Portico, an Open Source RTI HLA 1516e compliant. The hybrid HLA/FMU works as a time regulating and time constraining federate. On the Papyrus instance side, the time is managed by the DEScheduler of the Moka engine.

The interactions between the two technologies can be divided into two parts. The first part is when Papyrus instance must receive a message. The second part is when it must send a message. At each nanosecond of the simulation, a function is called by the time scheduler to check if nothing commits in the event buffer. We added a function which check modifications in the FMU file outputs. If any modification is noticed, a new event is triggered in order to treat this event. If the Papyrus simulation has to send a message, the Moka extension asks to DEScheduler the local time and get in return the simulation time in nanoseconds after a conversion to milliseconds. This time is given to the FMU with the doStep() function.

Time management between Papyrus instances

The time synchronization between Papyrus instances is an important point. When an instance of Papyrus is locked by another one, its local clock must continue to fly. This waiting time must be simulated by the Moka extension. The federate interactions goal is to calculate this duration. Figure 6 illustrates a simple example of this time synchronization mechanism. A more complex example was setup in order to test synchronization mechanisms. As a result, figure 7 is a BPMN diagram describing this simulation scenario. P1 starts and unlocks P2 & P3. Secondly, P1 must wait for P2 and P3 to finish their own simulation. Another interesting point is that P3 releases P1 after P2, even if P1 must first wait for P2. In that case, P1 stays naturally locked and wait for P2 message. Then, P1.5 is immediately unlocked (more explanation later in the paper with Figure 9). According to Figure 9, we can track execution steps of Figure 7 model. This example permits to test delay algorithm with a specific study case. Indeed, P1 must be firstly unlocked by P2.24 and secondly by P3.34. When P1 enters in P1.3 lock, it asks to P2 about its current local time: starting time = 2 and stores this information for later.

P3 process is faster than P2. As a result, Figure 9 diagram shows that P3 send a release message at t=7, but P1 is not yet locked by P3. So Federate 1 save this message and does not impact P1 simulation.

Later, at t=10, P2 releases P1.4 lock. At this time, the federate first unlocks P3.34 and determines its delay: d = release time (10) -starting time (2). Moka engine of Papyrus instance 1 can unlock its local time and set delay time to 8.

Then subsequently, the execution process goes in a second synchronized task P1.5, but it is immediately released: d = release time (7) -starting time (12) < 0 so no delay. Then, the distributed simulation ends since each Papyrus instance has reached its final task.

CONCLUSION

Papyrus community tends, in upcoming versions, to open the software to co-simulation with the implementation of the FMI standard. Nevertheless, as FMI does not include time management. To overcome this limit, we have presented in this work, a Moka federate component extended to communicate with an HLA RTI. This first result is pushing the platform capabilities to new time orchestrated features. It will allow developers to implement synchronization points between Papyrus instances.

The distributed simulation is currently at the early stage of design and implementation. So far, it can only share simple tokens as resources. Much remains to be done such as implementation of performance evaluation. For instance, requesting an FMU file and exchanging messages through a federation is time consuming compared to a standard HLA based DS compliant tool. Another interesting future work would be to implement more complex shared resources between Papyrus instances. To reach this goal, we plan to involve semaphore mechanisms. Finally, since another goal of the work presented in this paper is dealing with interoperability at large (e.g. databases, matlab files), we trust that HLA and FMI standards will facilitate connection with other heterogeneous tools or components.

Figure 1 :

 1 Figure 1 : Papyrus + Moka extension connected to FMU

Figure 2 :

 2 Figure 2 : Sequence flow diagram of Papyrus execution VS Papyrus execution with custom risks management profile

Figure 3 :

 3 Figure 3 : Custom UML Profile

Figure 4 :

 4 Figure 4 : HLA FOM file

Figure 5 :

 5 Figure 5 : Simple example of two Papyrus in BPMN

Figure 6

 6 Figure6explains the main exchanges between the components of the simulation. In our example, each Papyrus instance has two normal tasks (e.g. P1.1, P1.2 & P2.4, P2.5), and two synchronized tasks (P1.3 & P2.6). Those two last tasks extends a custom UML profile.

Figure 6 :

 6 Figure 6 : Sequence flow diagram of two Papyrus with synchronized point

Figure 7 :

 7 Figure 7 : BPMN process of 3 interconnected Papyrus

Figure 8 :

 8 Figure 8 : Logical architecture of Papyrus Distributed Simulation

Figure 9 :

 9 Figure 9 : Sequence flow diagram of Figure 5 BPMN Process



  Coordinated time, each federate advancing time is coordinated with each other  On messages, when the advance in time of a federate is triggered by reception of a message over HLA communication mechanism

 Optimistic, where each federate can advance in time at its own speed independently from the others.

Often used when the federate´s time is real time. If not, it receives a message from the past and triggers a rollback of the simulation.

In literature, authors worked on merging the HLA and the FMI standards. Yilmaz et al. presented in their paper (Yilmaz et al. 2014) a method to automatically generate an FMU entity able to join an HLA federation as member. A wrapper reads the FMU specification and generates the HLA layer. To demonstrate the FMUFd (FMU-Federate) usage, authors developed a simple DS example with MAK HLA RTI. In this paper, time synchronization between the two standards (FMI & HLA) is solved by updating the federate time at each running step of the FMU model.

AUTHOR BIOGRAPHIES

SIMON GORECKI is MSc and Ph.D. Student at University of Bordeaux in IMS Lab. His research domain is about simulating process with distributed simulations and HLA (High Level Architecture). His email address is simon.gorecki@u-bordeaux.fr JUDICAL RIBAULT is Ph.D. and part time Researcher at University of Bordeaux and IMS Lab. His research interests are methodologies for modeling and simulation.

GREGORY ZACHAREWICZ is Full Professor at IMT -Mines Ales and LG2IP Lab with both competences in enterprise engineering, discrete event modeling & distributed simulation. His email address is gregory.zacharewicz@u-bordeaux.fr NICOLAS PERRY is Full Professor at ParisTech ENSAM of Bordeaux. Hi domain of research is about system engineering, product process integration and green manufacturing. His email address is nicolas.perry@ensam.eu YVES DUCQ is Full Professor at University of Bordeaux in IML Lab. His domain of research is about enterprise modeling and process performances .His email address is yves.ducq@u-bordeaux.fr


E Duran Venegas

S Le Dizes

Generalised helical vortex pairs

New solutions describing the interaction of helical pairs of counter-rotating vortices are obtained using a vortex filament approach. The vortices are assumed to have a small core size allowing the calculation of the self-induced velocities from Biot-Savart law using the cut-off theory.

These new vortex structures do not possess any helical symmetry but they exhibit a spatial periodicity (figure 1) and are stationary in a rotating and translating frame.

Their properties, such as radial deformation, frame velocity (figure 2) or induced flow, are provided as a function of the four geometric parameters characterizing each solution:

Approximate solutions are also obtained when the mutual interaction is weak. This allows us to provide explicit expressions for the rotation and translation velocities of the structure in this limit. The variation of the vortex core size induced by the helix deformation is also analysed. We show that these variations have a weak effect on the shape and characteristics of the solutions, for the range of parameters that we have considered.

R ⇤ = R int R ext , h ⇤ = h ext R ext , ↵ = h int h ext , " = a Rext .
The results are finally applied to rotor wakes. It is explained how these solutions could possibly describe the far wake of an helicopter rotor in vertical flight (figure 3). velocity fluctuation with h ⇤ and R ⇤ , in agreement with the increase of the vortex deformation (see figure 10).

The mean tangential velocity of each vortex is shown in the (↵, h) plane in figure 15. A positive value corresponds to an advection in the positive axial direction, a negative value to an advection in the opposite direction. Not surprisingly, the mean tangential velocity blows up as ↵ ! 1 like ⌦ F and W F . It is also interesting to note that the contour curves are similar (in shape) for both vortices and close to those obtained for ⌦ F R 2 ext / in figure 13(a).

In figure 16a, we have displayed on the same plot the parameters for which mean tangential velocities and ⌦ F vanish. We clearly see that mean tangential velocities and velocity fluctuation with h ⇤ and R ⇤ , in agreement with the increase of the vortex deformation (see figure 10).

The mean tangential velocity of each vortex is shown in the (↵, h) plane in figure 15. A positive value corresponds to an advection in the positive axial direction, a negative value to an advection in the opposite direction. Not surprisingly, the mean tangential velocity blows up as ↵ ! 1 like ⌦ F and W F . It is also interesting to note that the contour curves are similar (in shape) for both vortices and close to those obtained for ⌦ F R 2 ext / in figure 13(a).

In figure 16a, we have displayed on the same plot the parameters for which mean tangential velocities and ⌦ F vanish. We clearly see that mean tangential velocities and

Figure 1 :

 1 Figure 1: Parameters defining the generalized helical vortex pairs.

Figure 2 :

 2 Figure 2: Angular velocity (right) and axial velocity (right) of the frame where the solution is stationary as a function of α and h * for R * =0.5 and ε=0.03. Both quantities have been normalized using the external radius Rext , and the vortex circulation Γ.

Figure 3 :Figure 13 :Figure 14 :

 31314 Figure 3: Diagram of the different rotor flow regimes. White regions: windmill brake regime or wind turbine regime (both vortices are going upwards as the external wind).Light gray regions: ascending regime (both vortices are going downwards as the external wind). Dark gray regions: slow descending regime and vortex ring state (both vortices are going downwards while the external wind is going upwards). Close to the line F = 0, there exists a small region where one vortex is going upwards while the other is going downwards: such a solution cannot describe the (far) wake of a rotor.

Figure 13 :Figure 14 :

 1314 Figure 13: Contour values of the frame velocity in the (↵, h) plane for R ⇤ = 0.5, " = 0.03, N = 1 and  = 1. (a) Angular velocity ⌦ = ⌦ F R 2 ext / , (b) Axial velocity W = W F R ext / . The dashed line (↵ = 1) indicates a line where ⌦ and W are not defined. The thick solid curve corresponds to the level zero.