François Queyroi

Comparing Static and Dynamic Graphs built from Mobility Traces A case study with maritime transportation data

Keywords: sequential data analysis, dynamic graphs, time-dependent shortestpaths, transportation networks

Introduction. This paper introduces ongoing researches on graphs built from aggregation of sequential data. More precisely, we focus on spatially and temporally located sequential data such as the records of ports visited by transportation vessels. The traces used here are extracted from the Lloyd's List Intelligence database 1and consists of container ships movements between April 1st 2009 and July 31th 2009 (made available by the WorldSeastems ERC project [START_REF] Ducruet | ERC World Seastems project[END_REF][START_REF] Ducruet | Mapping the globe. the patterns of mega-ships[END_REF]): the ports visited along with the arrival and sailing dates of the vessels. Although dynamic in nature, the shipping network is often studied as a static graph built by aggregation over a given period [START_REF] Hu | Empirical analysis of the worldwide maritime transportation network[END_REF][START_REF] Wang | A study of the temporal robustness of the growing global container-shipping network[END_REF]. For instance, [START_REF] Wang | A study of the temporal robustness of the growing global container-shipping network[END_REF] compared two static aggregations of transportation companies lines in 2004 and 2014. Looking at the robustness of some graph measures, they concluded that the maritime network became more vulnerable to port failures. However, there are several ways to construct a network from sequential data. They actually correspond to different aggregation strategies which lead to different interpretations of a same graph measure. The work presented here focuses on the comparison between network statistics depending on the chosen strategies. One important contribution of this work is the comparison between static shortest-path lengths and time-dependent fastest-path durations.

Aggregation strategies. In the context of maritime transportation networks, we use the concepts of space L and space A (using a terminology close to [START_REF] Hu | Empirical analysis of the worldwide maritime transportation network[END_REF]). For a vessel with a sequence of visited ports (a → b → c), a space L network contains direct links (a → b) and (b → c) while a space A network also contains the indirect link François Queyroi LS2N, UMR CNRS 6004, Université de Nantes E-mail: francois.queyroi@univ-nantes.

(a → c). The out-degree of a vertex (port) a in space L (respectively space A) corresponds to the number of direct destinations (resp. reachable ports) using a ship leaving a. The unweighted shortest-path distance between ports a and b in space L (resp. space A) corresponds to the minimum of number of visited ports (resp. number of ships used) from a to b.

We can also use knowledge about the domain to define networks in space P. Indeed, vessels are not independent agents but are assigned to a set of ports (lines or routes) by shipping companies. The assignations may change over time. Other researches [START_REF] Wang | A study of the temporal robustness of the growing global container-shipping network[END_REF][START_REF] Hu | Empirical analysis of the worldwide maritime transportation network[END_REF] directly use the lines given by shipping companies to define space P. For a given line, space P includes the indirect links between all ports visited on this route. Unlike space A, it does not include indirect links between ports that are never part of the same line and where no exchange of goods could have taken place. In space P, the out-degree corresponds to the number of lines available in a port while distances correspond to the minimum number of different lines required to go from a port to another. However, the Lloyd's database does not include shipping companies' planned lines but "only" the temporal sequences of ports actually visited by each vessel. The sequences still contain some patterns that can be used to recover the lines. For instance, a classic shipping route passes through a given port at most twice [START_REF] Shintani | The container shipping network design problem with empty container repositioning[END_REF] (following a circle pattern). We develop a simple procedure to extract lines from the sequence of ports visited by a given ship: we recursively cut the port sequence using as pivot the most visited port until the sequence is partitioned into valid circle sub-sequences. Notice that temporal and sequential information are lost using space L, P or A. However, we can use ships' arrival and sailing dates in each port and set the average time needed to cross an edge (a → b) as the average difference between sailing dates from a and arrivals dates in b. But this definition of edges weight still discards an important feature of temporal networks namely the frequencies of interaction between elements (number of ships going between two ports on a given period).

Another aggregation that preserves frequencies is the use of dynamic graphs model. The idea is to store every arrival and sailing dates for each edge in space L,P or A. The notion of shortest-path distance can be intuitively interpreted the fastest-path duration which can be computed using a variation of Dijkstra's algorithm [START_REF] Dean | Shortest paths in fifo time-dependent networks: Theory and algorithms[END_REF]. Notice that a path is here defined not only by a starting port but also by a starting date. This aggregation a priori also discards sequential information as one path can jump from one ship to another without any cost. This hidden hypothesis leads to strange results as a fastest path may include a lot of ports. We can introduce a stopover cost using multi-modal fastest path durations [START_REF] Delling | Accelerating multi-modal route planning by accessnodes[END_REF]. One solution here is computing fastest-paths using sailing and arrivals dates on edges in space P or A adding a fixed time c (time required to take another line or ship) each time a port is visited. Indeed, in space P (resp. A), if a fastest path contains an internal node it means that a change of line (resp. ship) had to take place. The additivity of shortest-paths is therefore not violated and Dijkstra's algorithm variation can still be used.

Results

We compare the four networks and computed distances using as vertex statistic the average distances from one port to the others. A low value indicates that a port is easily reachable. As the dynamic network may not be "temporally connected" we filtered out ports with a low frequency of interactions and took the average of travel duration over the period. We therefore compare the results on 438 ports. The density of connections for space L,P and A are 0.03, 0.1 and 0.15 respectively. First results suggest that space P and space A are actually very similar w.r.t. static distances and almost identical w.r.t. temporal distances. This may indicate that most ships follow regular routes on this relatively short period. They also suggest that unweighted distances in space L (i.e. minimum number of visited ports) is a better approximation for travel durations than the unweighted distances in space P or A (minimum number of used lines or ships). Although correlations with space P and A distances exist, it indicates that assessing ports accessibility only using space P or A may discard important information. This effect is stronger when looking at average travel durations with a cost of stopover c = 2. These correlations seem weaker when looking at the weighted distances (i.e. using the average time needed to cross an edge). This observation is counter-intuitive and further tests are needed to validate it.

Ongoing work and future directions. More work is needed in order to assess the effects of errors in the database. Indeed, some trips between distant ports are shorter than they should. This may affect weighted distances and fastest-paths computation. Moreover, our definition of space P relies on one possible route extraction algorithm that should be compared against others or directly against the shipping lines given by the shipping companies. Future directions include the analysis of the measurements and their robustness to change on different time periods [START_REF] Wang | A study of the temporal robustness of the growing global container-shipping network[END_REF]. We will also compare other network measures (e.g. clustering coefficient) with their dynamic counterparts (e.g. a dynamic clustering coefficient corresponding to a local and dynamic version of betweeneess centrality).

Data and algorithms' implementations used are available at https://github.com/ fqueyroi/tulip_plugins/tree/master/papers/NetFromSequentialData


Julen Munárriz

Rubén Laplaza

A Martín Pendás

Julia Contreras-García

A first step towards Quantum Energy Potentials of Electron Pairs

A first step towards the construction of a quantum force field for electron pairs in direct space is taken. Making use of topological tools (Interacting Quantum Atoms and the Electron Localization Function), we have analyzed the dependency of electron pairs electrostatic, kinetic and exchange-correlation energies upon bond stretching. Simple correlations were found, and can be explained with elementary models such as the homogeneous electron gas. The resulting energy model is applicable to various bonding regimes: from homopolar to highly polarized and even to non-conventional bonds. Overall, this is a fresh approach for developing real space-based force fields including an exchange-correlation term. It provides the relative weight of each of the contributions, showing that, in common Lewis structures, the exchange correlation contribution between electron pairs is negligible. However, our results reveal that classical approximations progressively fail for delocalized electrons, including lone pairs. This theoretical framework justifies the success of the classic Bond Charge Model (BCM) approach in solid state systems and sets the basis of its limits. Finally, this approach opens the door towards the development of quantitative rigorous energy models based on the ELF topology.

Introduction

The simulation of big systems requires the use of analytical force fields (FFs) that enable the performance of molecular dynamic simulations in reasonable computing times. This is usually done by following simple atomic models. For example, bond lengths are commonly analysed in terms of the Hooke's law, where potential and kinetic energy balances are parametrized to adjust to experimental/computational results. 1,2 Nonetheless, in many cases, chemical behaviour is determined by non-classical terms (e.g. exchange or correlation), which are not easily encoded in these potentials. We propose hereby a change of paradigm that is able to easily capture electron reorganization (towards building reactive FFs), as well as quantum terms.

Atomic models usually require specific parametrization for different coordination numbers. This stems from the fact that electrons (i.e. bonds) are the ones being re-allocated during a reaction. Therefore, the static view of electrons belonging to one atom, and hence, its oxidation state, can be avoided if the electron pair is taken as the main unit. It should be noted that this approach has been used for a long time in chemistry, although not in a quantitative manner. Before the advent of formal force fields, the use of mere electrostatic reasoning among electron pairs allowed for the first time to predict molecular structures, in what is nowadays known as the Valence Shell Electron Pair Repulsion (VSEPR) theory. 3 Moreover the use of electron pairs also provides a simpler link to quantum mechanics. The usefulness of electron pairs as the basis of the representation lies on the fact that paired electrons behave very much like bosons. 4,5 This means that Pauli correlation is minimized in a system composed by perfectly localized electron pairs, so that the exchangecorrelation term (more difficult to disentangle with simple models) is minimal. These advantages have been exploited in both, classical approaches, such as the Bond Charge Model (BCM), 6,7 where only potential and kinetic energy terms are included; as well as in more updated formulations based on quantum mechanical calculations. 8 In order to provide robust and rigorous analytical expressions for a quantum FF, it is necessary to count with a recipe for the identification of electron pairs, and their parametrization. For that, the Interacting Quantum Atoms (IQA) energy decomposition scheme is very suitable, since it allows to obtain and quantify the different energetic contributions from the real space partitions of the system. [START_REF] Pendás | Two-electron integrations in the quantum theory of atoms in molecules[END_REF]10 The advantages of the IQA scheme are showcased by its success in the study of different chemical phenomena, such as halogen and hydrogen bonding, 11,12 and metal-ligand interactions, 13 among many others. However, these calculations are expensive and limited to small systems. Therefore, by reconstructing models from this scheme, the advantages of IQA could be transferred to bigger systems.

For translating the IQA partition to electron pair energy models, an electron-pair real-space partition is needed. The IQA scheme is usually coupled to the Quantum Theory of Atoms in Molecules (QTAIM) partition in order to associate the different energy terms to atomic contributions. In this sense, it is worth 14 This framework also proved to be suitable for geometry optimizations when combined with machine learning approaches. 15,16 However, since we want to focus on electron pairs, we will couple IQA to the Electron Localization Function (ELF) topology. 17 The ELF is able to recover Lewis entities, as bonds and lone pairs, providing the classical VSEPR picture, and is thus the perfect starting point for our energy model, allowing to quantify (classic and quantum) interactions between electron pairs. 18 In this article we propose a first step in constructing a stretching energy model based on the ELF topology. In order to do so, we evaluate electron pair interactions through classical (kinetic and electrostatic) and quantum energy terms (exchange and correlation). The ease and chemical rigor with which the model parameters can be obtained set the foundations of a fresh approach to real space force fields. In addition, a careful analysis of the relative weight of the classical and quantum energy terms is presented. Valuable implications in general force field development are also derived: we identify the chemical bonds that require quantum terms among those that can be studied under a purely classical perspective. The resulting model is beautifully connected with previous "intuitive" models (e.g. VSEPR, BCM). It also sets the basis for connecting the ELF topology with the energy of the system, leading to a simple chemical model that attempts to answer the initial question: "how do electron pairs (bonds, lone pairs) contribute to the energy of the system?"

Theoretical and Computational Methods

All energy scans and wavefunctions for ELF and IQA analyses were obtained through DFT calculations, using the Gaussian 09 program package. [START_REF] Frisch | Gaussian 09, Revision D.01[END_REF] The B3LYP exchange-correlation functional [START_REF] Lee | [END_REF]21 was used, in conjunction with the 6-31G(d) basis set. 22,23 Stability checks were carried out with other functionals (PBE and PBE0) [24][25][26] and with Hartree-Fock, in conjunction with bigger basis sets (6-311G(d,p), cc-PVDZ and cc-PVTZ). 27 The results obtained by varying the computational conditions led to the same behaviour, yielding very similar fittings, while the computational cost increased considerably with the increase in the basis set size. These results are provided in the Supporting Information (see Tables S1-S6). Hence, the small basis set was favoured. Bonding descriptors were computed at the B3LYP/6-311G(d,p) level to further ensure basis set convergence of the density at the bond critical points.

The partition into localized electron pairs was carried out with the Electron Localization Function (ELF). 17 The ELF kernel, , can be interpreted as a measure of the excess of local kinetic energy due to the Pauli Principle, relative to the homogeneous electron gas kinetic energy density. 28 The ELF,  , is mapped through a Lorentzian function (see eq. 1), to a scale ranging from 0 (when  → ∞), to 1 (when  → 0).

(1)

The gradient of this function, , is used to induce a topological partition which divides the space into non-overlapping regions (basins) whose properties can be determined by integrating the appropriate densities over their associated volume. Hence, if one is interested in, for example, lone pairs populations, it suffices to integrate the electron density, , over the corresponding region associated to the lone pair maximum.

The ELF studies were performed with a locally modified version of the TopMod program, 29 using the monodeterminantal B3LYP/6-31G(d) wavefunctions in conjunction with a cubic tridimensional grid of 200 points in each direction. ELF plots were made with Chimera software. 30,31 The energies of the different topological basins can be calculated with the Interacting Quantum Atoms (IQA) energy decomposition scheme. 32 This approach provides a set of unique and rigorous energetic terms that additively recover the exact energy of the system. Unlike many topological analyses, this method is not only suitable for stationary points (e.g. equilibrium geometries), such as with virial related energy partitions, but also for non-equilibrium geometries. This feature was crucial for evaluating the energy terms along bond elongations in geometry scans. The energy terms are calculated by partitioning the first and second-order density matrices with respect to real space partitions, usually the QTAIM atomic basins (see eq. 2).

∑ ∑ ∑ ∑ (2)
This way, the total energy is divided into an intra-atomic contribution,

, and the inter-atomic interaction between each pair of atoms, . The former is the sum of the kinetic energy of electrons in atom A, , and the electron-electron and electron-nucleus interactions, and , respectively. is the sum of the interaction among the electrons in A and the nucleus of B, ; the nucleus of A and the electrons in B,

; and the nucleus-nucleus and electron-electron interactions between A and B, and , respectively. The Electron-electron repulsion can be further divided into a classical (electrostatic) contribution, , and a non-classical one, , which is the sum of exchange (X) and correlation (Corr) terms (see eq. 3).

(3) Following this scheme, each contribution to the total energy (intra and interatomic) can be expressed just like in atomistic approaches, that is, as a sum of kinetic and electrostatic terms, plus a non-classical interatomic term, as shown in eq. 4 and 5.

(4) (5) Where all electrostatic terms have been put under a common "Coul" index: and . IQA can also be applied to ELF partitions. When considering bonding basins (and valence in general) within the IQA partition, the nuclear terms presented in the previous equations become null ( 0, 0, 0). Therefore, for a bonding basin, "Bond", which interacts with a core basin that represents an atom A, the different energy terms can be expressed as shown in eq. 6 and 7.

(

Since we are interested in developing an energy model accounting for interactions between electron pairs, this was the selected approach. This way, it is possible to calculate interactions between atoms and "bonds", as well as interactions between two "bonds" and "bonds" with "lone pairs". For this purpose, the original code was modified so as to perform integration tasks over ELF basins. 18 The IQA-ELF approach provides an accurate reference that can be used to analyse the behaviour of the energy terms and to construct energy potentials that take into account classical and non-classical terms. In the next sections, we will examine the dependency of these terms with respect to bond elongation and compare them with previous potentials.

The original IQA implementation could only deal with HF wavefunctions; however, recent developments provide support for DFT-derived ones. 10 IQA calculations were performed with PROMOLDEN. [START_REF] Pendás | Two-electron integrations in the quantum theory of atoms in molecules[END_REF]33 Results and Discussion

The model

As a starting point, the ELF-IQA approach was applied to a series of simple yet representative molecules in the CH 3 -X (X = CH 3 , NH 2 , OH, F) series. In order to analyse the variation of the different terms that contribute to and (A-Bond), and attempt to rationalize the observed behaviour; the C-X interatomic distance, R, was stretched around the equilibrium position, R eq , approximately in the range [0.90-1.10]•R eq . Since we are resorting to small elongations, the use of a DFT approach is justified. The following standard ELF notation is used: V(C,X) represents the C-X bond, V(X) the lone pairs on X, C(Y) or simply Y, designates the core of atom Y and so on. For the sake of clarity, we will only analyse the results obtained for CH 3 -NH 2 in detail, as an illustrative example with bonds and lone pairs. Similar results have been obtained for all the molecules in the CH 3 -X test set, and are provided in the Supporting Information.

The CH 3 -NH 2 ELF basins and their populations are shown in Figure 1a. The evolution of the ELF populations of V(N) and V(C,N) upon bond stretching follows the expected behaviour, with a flux of electrons towards the lone pairs as the distance increases (see Figure 1b). In order to develop working equations along the stretching, we need to understand the evolution of the main parameters during this process. Bond charges and bonding basin sizes were found to change linearly with R (as shown in Figures 1b and1c). Similar results for the complete CH 3 -X set are provided in Figures S1 andS2.

In the following sections, the interacting potentials for the kinetic, electrostatic and exchange correlation terms involving the C-X bond and valence lone pairs will be developed from the ELF-IQA partition. In all cases, the difference between the calculated IQA energy and the DFT obtained one (∆E) was used to evaluate the numerical integration accuracy. The error held roughly constant and lower than 1% of the total energy, as shown in Table S7.

Classical electrostatic interaction energy. Since the ELF basins are non-overlapping, the electrostatic interaction between them can be reasonably well approximated by the Coulombic interaction between point charges located at the barycenter of the basin charge density. We have carried out a test of the performance of such an approximation in the reference molecule, CH 3 -NH 2 . The centre of charges for the C-N bond only deviates 0.07 a.u. from the ELF attractor. This small difference is clearly negligible when compared to the bond length, 2.75 a.u. Hence, monopole electrostatics are expected to reproduce reasonably well the electrostatic interactions (see eq. 8).

∝ (8)

Where q i refers to the electronic charge of basin i (be it core, bond or lone pair). Hence, the overall electrostatic interaction should follow the classical interacting potential Eelec = B/R + C, where B and C are constants. For testing this hypothesis, we used the ELF-IQA approach to calculate the classical component of and (see eq. 6 and 7) for the interactions between the C-X basin (V(C,X)) and both the C and X (C(C) and C(X)) ones. In addition, we computed the interactions between the X lone pairs (V(X)) and X (C(X)). As previously explained, these interaction between basins represents the interaction between the C-X bond and C and X atoms, as well as the interactions between the lone pair in X and the X core, respectively. We represented these contributions against 1/R. Results for the interaction between the pairs of valence basins in methanamine are shown in Figure 2 (see Figure S4 for the results involving all other molecules). As it can be seen, the linear regression fittings present very good correlation coefficients, with R 2 being higher than 0.97 in all cases, for both intra and inter-basin contributions. The good agreement confirms that the ELF maximum can be used to approximate a point charge behavior. Remarkably, the linear correlation between the electrostatic interactions and 1/R also holds for intra-basin interactions, where q A = q B (i.e. self-interaction can be neglected in the energy-model). Note that repulsive interactions are obtained for intra-basin terms, since they correspond to electron-electron repulsions, whereas interbasin terms are positive, due to the extra attractive contribution from the nucleus. This agreement is also maintained when the basin charges (i. e. the electron density integrated over the basin volume) are taken into account (see Figure S5). Overall, this indicates that zero-th order electrostatics is a reasonably accurate approach for ELF basins. Kinetic energy. In order to analyse the evolution of the valence electron pairs kinetic energy, we can assume that the electron density in the low density region (i.e. bonds far from the core) can be described by a zero-th order expansion of the kinetic energy density of a homogeneous electron gas, / .

Assuming now a linear behavior of charges upon stretching (as shown in Figures 1b andS2), we obtain a simple expression relating the kinetic energy of valence basins with distance (see eq. 9). The development of eq. 9 is provided in the Supporting Information.

/

The expression provided in eq. 9 depends on the bond basin radius (R B ). In order to verify whether we can relate it with the C-X bond distance, we calculate R B as a function of R. For that, we subtracted the core radii to the interatomic distance, as proposed by Savin et al. 34 The results show that R B (and R LP ) are proportional to R (see Figure S3). Therefore, we can develop a working expression of the form: T ∝ R -2 . The validity of this relation for the kinetic energy of the valence basins of CH 3 -NH 2 , V(C,N) and V(N), is proven in Figure 3. The fittings for other molecules are provided in the Supporting Information (Figure S6). Overall, the agreement within the whole series of molecules leads to regression coefficients, R 2 , ranging from 0.984 to 1.00 for both, bonds and lone pairs. Exchange-Correlation energy. Finally, let us examine the quantum mechanical contributions. Exchange-Correlation interactions have been shown to follow a linear behavior with the delocalization index, δ AB , according to eq. 10. 35 (10) Delocalization indices between QTAIM atoms have been shown to decay exponentially with the inter-centre distance for insulators, 36 as they depend on orbital overlap. In the case of ELF partitions, we have found that δ AB has a clear hyperbolic behavior (that is, δ AB = k/R) for interactions between the bond and core basins in the range of distances considered (as shown in Figure S7). Accordingly, the effective expression for the exchange-correlation energy should follow eq. 11. Figure 4 shows how well this relation holds for the amine, the fittings for all the other molecules being provided in Figure S8 of the Supporting Information. Very good linear regression coefficients, R 2 > 0.99 in most cases, were obtained, confirming the validity of the previous reasoning and the derived dependencies. (11) Interplay of terms. The previous results have important consequences that are highly relevant in both, force field development, and quantum chemical topology. Firstly, they confirm the ability of an electron pair force field (that includes chemical bonds and lone pairs) derived from topology to describe a system while using well-behaved quantum contributions. The latter are easy to parametrize and to be added to available "classical" potentials.

Secondly, the quantum mechanical exchange-correlation energy term is found to be noticeably smaller than the classical terms for all systems in the dataset (see Figure 5). As an example, we provide the different contributions for methanamine in the equilibrium position. For the V(C,N) basins:

T V(C,N) = 2.207 a.u., E elec-intra V(C,N) = 0.825 a.u., E elec-inter V(C,N)-C = -4.028 a.u., E elec-inter V(C,N)-N = -5.723 a.u., E x-c V(C,N)-C = -0.044 a.u., E x-c V(C,N)-N = -0.058 a.u. For V(N): T V(N) = 3.281 a.u., E elec- intra V(N) = 1.401 a.u., E elec-inter V(N)-N = -9.379 a.u., E x-c V(N)-N = - 0.138 a.u.
It is interesting to note that exchange-correlation terms are two orders of magnitude smaller than classical terms for bond basins. Nonetheless, they are only one order of magnitude smaller in the case of lone pairs. This reveals the greater relevance of quantum effects in the description of lone pairs. When building an electron pair force field they can be neglected, in a first approximation, in the case of covalent bonds, but should be taken into account when developing accurate interacting potentials for lone pairs. Furthermore, this prevalence of quantum effects further justifies the explicit treatment of lone pairs in Force Fields, which are often simply included in atomic contributions.

As a final check, Table S8 shows that intra-basin core energies (kinetic, electrostatic and exchange-correlation) remain constant upon elongation. Hence, they can be ignored in the construction of the interacting potential.

Comparison with previous models

Up to now, we have coupled IQA with the ELF topology to obtain the different interaction energy terms involving Lewis-like entities (atoms, bonds and lone pairs). According to the previous discussion, the total bond energy would be the sum of kinetic, electrostatic and exchange-correlation contributions, even though in most cases the latter will play a minor role. Interestingly, this framework is formally equivalent to the one provided by the Bond Charge Model (BCM). 6 This approach is a simple albeit chemically intuitive model for calculating the energy of molecules (or solids). In particular, a homonuclear diatomic molecule, say A 2 , is represented as two cores with +q/2 positive charge, each of them interacting through a bond holding a charge of -q electrons (see Figure 6 top). The bond charge is allowed to move along the bond length (R B ), which is a fraction of the equilibrium interatomic distance (R), R B = υR. The total energy, which depends on the charge and the bond length, is the sum of three different contributions (see eq. 12, in atomic units).

, 2 D′ (12)
Where E 0 is the core energy (equal to 2E A for a homonuclear molecule, with E A being the core energy of each atom); E 1 accounts for the classical electrostatic interactions (core-bond and core-core); and E 2 represents the kinetic energy of the bond electrons moving along the bond length, R B . The model can also be expanded to heteronuclear systems, say AB, (see Figure 6 bottom) leading to eq. 13. This expression is formally equivalent to the one obtained for homonuclear molecules.

, ′ D′ (13)
Originally, the model parameters were obtained by fitting to experimental data. 6,7 However, it is remarkable that the BCM model has the same energy terms that have been analysed by means of the ELF-IQA approach (with the exception of the exchange-correlation term, which is non-classical) and that they present the same dependencies with the interatomic distance, R. The parallelism with the terms derived from our ELF-IQA energy model is thus evident:

′ ′ 0
The first three equations have been shown to work for our test set. As for the exchange-correlation contribution, a comparison of the scales of Figures 3 and4 provides a first idea of the error. As previously stated, bond exchange-correlation terms are typically two orders of magnitude smaller than the other terms for bonds and one order of magnitude for lone pairs (see Figure 5).

It is important to note that this framework explains the previous success of the BCM model, as the energy terms with a classical origin are the only ones that contribute significantly to the scaling law of the total energy of molecules whenever electrons are well-localized. It also justifies the success of ELFbased parametrizations of the model. Moreover, it explains why this model was successful for extended solids in material science, while not so much applied for open structures (layers, chains), where lone pairs are usually present. 37 As we have seen, neglecting the exchange-correlation terms is straight forward for bonds (two orders of magnitude smaller than kinetic and electrostatic terms), but not so much for lone pairs (only one order of magnitude smaller). The next section will be devoted to the analysis of the limits of the IQA-ELF potential.

Limits of the model

Since the model is based on localized electron pairs, it is perfectly compatible with the description of covalent bonds. Accordingly, it may not suffice to capture other types of interactions. Several examples pertaining to distinct bonding regimes have been analysed as an attempt to estimate the range of applicability of the model. On the basis of the extent of electron pair localization, it should be possible to predict outlier cases where the IQA-ELF potential should fail. The following systems have been addressed: BH 3 -NH 3 , BeH 2 -NH 3 , CH 3 -Li and Li 2 , as paradigmatic examples of dative covalent bond, noncovalent beryllium bond, ionic and metallic systems (see Figure 7). 38,39 It was found that both, dative and very polar bonds, still follow the proposed dependencies for the IQA-ELF potentials. Specifically, the BH 3 -NH 3 and BeH 2 -NH 3 molecules lead to linear regression coefficients, R 2 , for the electrostatic and bond kinetic energies higher than 0.97 and 0.99, respectively (see Table S9 andS10). The exchange-correlation contribution was also found to be negligible in comparison with the classical energy terms, being two orders of magnitude lower, as in the CH 3 -X series. Similar results are obtained for CH 3 -Li. Interestingly, the Li-V(C,Li) electrostatic term is smaller than the C-V(C,Li) counterpart by almost one order of magnitude. This result suggests that this contribution is not energetically meaningful, which can be related with the partially ionic character of the bond, which leads to the bonding basin being significantly closer to the carbon than to the lithium core (0.829 Å for C vs 1.169 Å for Li). The proposed scaling laws of the energy terms do not hold for Li 2 . From the visual point of view, a rather remarkable feature is already identifiable in Figure 7d. The volume of the bond basin (where a non-nuclear maximum of the density appears) is much bigger than in previous cases (898.23 vs 17.91 Å 3 in ethane). Moreover, the population in this basin also follows a different behavior upon stretching. It remains rather constant (between 1.81 and 1.83 e-) independently of the Li-Li interatomic distance. The resulting ELF-IQA kinetic, electrostatic and exchange-correlation energy terms are shown in Figure 8. It is noteworthy that all the terms are very small, of the order of 10 -4 a.u., very different from the ones obtained for CH 3 -X molecules, which were of the order of 10 1 a.u. (kinetic and electrostatic terms) and 10 -1 a.u. (exchange-correlation term). We can thus consider that the model works properly for all the molecules considered but Li 2 , for which it is not suitable.

Overall, the model is shown to be valid for covalent and partially covalent links, while failing for metallic bonds (as the case of Li 2 ). It should be noted that these two categories can be easily identified from properties derived from the electron density. Over the years, numerous attempts have been made to provide a distinct classification of the chemical bond from topological characteristics, 40,41 with special emphasis on the metallic bond. [42][43][44] Complex descriptors, as the metallicity index (ξ m ) proposed by Ayers et al. 44 have been used to further collect the character of bonds through topological measurements. This index is constructed to highlight regions of planar density, its expression being provided in eq. 14.

⁄ /

High values of ξ m are indicative of metallic character (typically ξ m > 25 in solids); while weak metallic bonds are identified by ξ m values between 1 and 5. When ξ m is inferior to 1, the bond presents a non-metallic character. Bonds of all the previous molecules were characterized through their electron density, ρ(r), its Laplacian, 2 ρ(r), and the aforementioned derivative descriptor at the bond critical point. Given that local topological indices extrapolate the characteristics of a single point to a large chemical entity, a larger basis set has been used to guarantee converged densities and robustness. Results are presented in Table 1.

First of all, it is clear that the limiting cases we chose lead indeed to positive Laplacians and larger values of ξ m . The metallicity index for Li 2 is much higher than for the rest of the considered bonds, which correlates with the breakdown of our model. Furthermore, the results do partially justify the lesser degree of accuracy of the model for the B-N, Be-N and C-Li bonds, as they exhibit greater ξ m indices and positive Laplacian values at their BCP. This feature can be related to their onlypartial covalent nature. All this is compatible with the basic idea encompassed here: localized pairs can be parametrized as point charges, so that when electrons indistinctly delocalize over a is not defined when 2 ρ < 0, yet it can be calculated to exemplify the significant difference between bonds. larger part of the system, the model fails. Contrarily, highly covalent bonds, with well localized pairs around their BCP can be assumed to be properly described by the model hereby presented.

Conclusions and outlook

We have presented a simple energy model based on the ELF topology for the stretching of bonds that explains the interaction between electron pairs while retaining their quantum mechanical character. Kinetic, electrostatic and exchange-correlation terms have been shown to follow distance dependencies that are akin to straightforward scaling laws. Electron pairs, as described by the Electron Localization Function, can be modelled from the homogeneous electron gas point of view. The proposed scaling behaviours have proven to work reasonably well for different polarities (from homonuclear to highly polar bonds, including non-conventional bonds). Whenever the electron pair fails to describe bonding, like in metallic bonds, or when highly delocalized bonding regimes set in, the model is expected to fail. We have also shown how simple local measures can be introduced to identify these situations.

Moreover, the energy model hereby introduced is able to explain the success of simple previous proposals, such as the Bond Charge Model, and its updated version, the ELF-BCM model. 37 In both of them, the quantum nature of the electrons was neglected. Our test set shows that, in general, the quantum mechanical contribution (represented by the exchangecorrelation term) is much smaller than the other terms for bonds. On the contrary, quantum contributions become more important as electron delocalization increases, as in the case of lone pairs. Provided that lone pairs are of significant importance in many chemical systems that are usually modelled with Force Field approaches, we hereby justify the careful consideration of their position and electrostatic behaviour. In this regard, the definition of explicit lone pairs has already been considered in molecular mechanics modelling, either mathematically or conceptually. 45,46 Notably, the development of a molecular energy model such as the one here devised may have important implications in other closely related fields, such as conceptual DFT (c-DFT). 47 The existence of a simple parametrization for the energy of a system would allow to obtain molecular properties (electronegativity, chemical hardness, etc.) in a simple manner and starting from chemically meaningful concepts. 48,49 As previously stated, another field on which this energy model may suppose a breakthrough is in Force Field development. In this respect, our model provides useful insight on non-classical terms and further justifies the design of Force Fields transcending the all-atom paradigm. Such Force Fields have already achieved remarkable success, 50 and we expect to foster these novel modelling methods. The work herein presented also opens the doors towards the further development of energy models based on Quantum Chemical Topology.

Figure 1 .

 1 Figure 1. a) CH 3 -X ELF basins (isovalue = 0.8) and populations (at the equilibrium geometry). V(C,H) basins are shown in blue, V(C,C) and V(C,N) in green, and V(N) in orange. b) V(C,N) and V(N) populations as a function of the C-N bond distance. Linear regression coefficients: R 2 = 0.982 for V(C,N) and 0.989 for V(N). c) V(C,N) and V(N) volumes as a function of the C-N bond distance. Linear regression coefficients: R 2 = 0.631 for V(C,N) and 0.988 for V(N).

Figure 2 .

 2 Figure 2. ELF-IQA computed classical electrostatic interaction energy (intra and inter) against 1/R for a) V(C,N) and, b) V(N) in CH 3 -NH 2 . Linear regression coefficients: R 2 = 0.978 for E intra V(C,N); 0.995 for E inter V(C,N)-C; 0.996 for E inter V(C,N)-N; 0.989 for E intra V(N); and 0.994 for E inter V(N)-N.

Figure 3 .

 3 Figure 3. ELF-IQA computed bond kinetic energy against 1/R 2 for V(C,N) and V(N) basins in CH 3 -NH 2 . Linear regression coefficients: R 2 = 0.997 for V(C,N) and 0.984 for V(N).

Figure 4 .

 4 Figure 4. ELF-IQA computed exchange-correlation energies against 1/R 2 for V(C,N) and V(N) basins with C and N in CH 3 -NH 2 . Linear regression coefficients: R 2 = 0.998 for V(C,N)-C; 0.997 for V(C,N)-N; and 0.988 for V(N)-N.

Figure 5 .

 5 Figure 5.Figure 5. ELF-IQA energy terms for C-N bond in CH 3 -NH 2 .

Figure 6 .

 6 Figure 6. Bond Charge Model representation of a homonuclear molecule (top) and a heteronuclear one (bottom).

Figure 7 .

 7 Figure 7. ELF basins (isovalue = 0.8) and populations for a) BH 3 -NH 3 , b) BeH 2 -NH 3 , c) CH 3 -Li, and d) Li 2 . Hydrogenoid, disynaptic and core basins are shown in blue, green and purple respectively.

Figure 8 .

 8 Figure 8. Kinetic, classical (intra and inter) and exchangecorrelation energy terms for Li-Li bond in Li 2 . For simplicity, energy values are multiplied by 10 4 .

  the previous work of Popelier et al., who developed the first protein FF based on the QTAIM topology.

	a. Departamento de Química Física and Instituto de Biocomputación y Física de
	Sistemas Complejos (BIFI), Universidad de Zaragoza, 50009, Zaragoza, Spain
	b. Sorbonne Université, CNRS, Laboratoire de Chimie Théorique, LCT, F. 75005 Paris,
	France.
	c. Laboratorio de Química Física y Analítica, Universidad de Oviedo, 33006, Oviedo,
	Spain.
	†Electronic	Supplementary	Information	(ESI)	available.	See
	DOI: 10.1039/x0xx00000x

noting

Table 1 .

 1 Topological characteristics of the selected bonds. Calculations were performed at the B3LYP/6-311G(d,p) level of theory.

	C-C	0.238	-0.532	-5.518
	C-N	0.261	-0.670	-5.116
	C-O	0.256	-0.504	-6.560
	C-F	0.233	-0.022	-130.6
	B-N	0.099	0.420	1.631
	Be-N	0.057	0.337	0.801
	C-Li	0.044	0.213	0.821
	Li-Li	0.013	0.002	12.78
	* ξ m

Bond

 BCP (a.u)

2

 BCP (a.u) ξ m ,BCP*

Acknowledgements

J.M. acknowledges financial support provided by the Spanish Ministerio de Ciencia, Innovación y Universidades (FPU14/06003 and EST17/00161), as well as from the Universidad de Zaragoza, Fundación Bancaria Ibercaja and Fundación CAI (CB 6/17). R.L. acknowledges ED388 for a PhD. grant. A.M.P expresses its gratefulness for the financial support provided by the Spanish Ministerio de Ciencia, Innovación y Universidades (CTQ2015-65790-P). J.C.-G. acknowledges CALSIMLAB and the ANR within the Investissements d'Avenir program under reference ANR-11-IDEX-0004-02. In addition, the resources from the supercomputer "memento", as well as the technical expertise and assistance provided by BIFI-ZCAM (Universidad de Zaragoza) and the Laboratoire de Chimie Théorique (Sorbonne Université) are gratefully appreciated.

Conflicts of interest

There are no conflicts to declare.


Mirzaie Majid
email: majid.mirzaie@znu.ac.ir

Ataabadi

Ali Bahrami
email: a.bahrami@sci.ui.ac.ir

Mehdi Yazdi
email: m.yazdi@sci.ui.ac.ir

André Nel
email: anel@mnhn.fr

A locust witness of a trans-oceanic Oligocene migration between Arabia and Iran (Orthoptera: Acrididae)

Keywords: Insecta, Acrididae, Formicidae, putative migration, Oligocene, Parathetys, Iran

Fossil insects are very rarely found in sediments of deep marine origin. Nevertheless they can be of great interest to trace past events such as trans-oceanic migrations. Here we document the first fossil insects from Iran, viz. several alate ants dead during mating swarms and a migratory locust, found in the pelagic Oligocene sediments of the Zagros Mountains. This locust represents the first putative indication of insect migrations between the Arabian-African and Asiatic continents through the Parathetys, probably in relation with the development of open grassland biotas in these areas.

Introduction

Direct fossil evidence of a long-distance mass insect migration are extremely rare. One record is known for a moth migration in the diatomous sediments of the Fur Formation from the Early Paleogene of Denmark (55 million years old) [START_REF] Rust | Fossil record of mass moth migration[END_REF]. The main reason of the rarity of such evidences is that insects need to be fossilized in marine sediments. The discovery of numerous fossils of the same species in the same layer of a terrestrial deposit can correspond to a continental migration but also to a local mass occurrence. For instance, several griffenflies were found in one layer of the Middle Permian of Lodève (France) but it is not possible to say that they correspond to a migration [START_REF] Nel | The Odonatoptera of the Late Permian Lodève Basin (Insecta)[END_REF]. [START_REF] Arillo | The fossil Acrididae from the Oligocene of Izarra (Alava, Spain). The antiquity of gregarious behavior (Orthoptera: Caelifera)[END_REF], on the basis of the numerous co-occurrences of locusts fossilized in the same Spanish Oligocene paleolake, hypothesized that they could correspond to an evidence of gregarious behavior, but they could not infer a phenomenon of mass migration. Insect taphocoenosis of marine origin are not very frequent. Notable exceptions are the Early Triassic of Australia [START_REF] Haig | Early Triassic (early Olenekian) life in the interior of East Gondwana: mixed marine-terrestrial biota from the Kockatea Shale, Western Australia[END_REF], Early Jurassic of Northern Europe [START_REF] Ansorge | Insekten aus dem Oberen Lias von Grimmen (Vorpommern, Norddeutschland)[END_REF], Late Jurassic of Germany [START_REF] Viohl | The paleoenvironment of the Late Jurassic fishes from the Southern Franconian Alb[END_REF] and France [START_REF] Nel | Exceptionally preserved insect fossils in the Late Jurassic lagoon of Orbagnoux (Rhône Valley, France)[END_REF], Early Cretaceous of Lebanon [START_REF] Prokop | First record of a fossil beetle (Coleoptera, Haliplidae) from the base of Paleocene flysch sediment in the Magura unit (Outer Western Carpathians, Moravia)[END_REF][START_REF] Vršanský | Mieroblattina pacis gen. et sp. n. -Upper Cretaceous cockroach (Blattida: Mesoblattinidae) from Nammoura limestone of Lebanon[END_REF], and Early Paleogene of D enmark [START_REF] Rust | Biostratinomie von Insekten aus der Fur-Formation von Dänemark (Moler, oberes Paleozän/unteres Eozän) [Biostratinomy of insects from the Fur Formation of Denmark (Moler, Upper Paleocene / Lower Eocene)[END_REF]. These relatively rich outcrops correspond to relatively shallow water, frequently lagoonal, paleoenvironments close to land masses. Sediments of pelagic origin normally do not contain fossil insects, with the exceptions of a beetle elytron found in Paleocene flysch sediments produced in the lower bathyal to abyssal zone of Moravia [START_REF] Prokop | First record of a fossil beetle (Coleoptera, Haliplidae) from the base of Paleocene flysch sediment in the Magura unit (Outer Western Carpathians, Moravia)[END_REF], and the recent discovery of a fly found in the deep-water marine Early Miocene of Zillerleite, Germany (Nel, pers. obs.).

Here we record the first fossil insects from Iran, viz. remnants of ants and of the wings of a locust, found by one of us (MMA) in the pelagic Oligocene sediments of the Zagros Mountains. This locust probably died during a migration between the Arabian Peninsula and the continental Iran. Although fossil insects were known from these fossil beds [START_REF] Garassino | Report on decapod crustaceans from the Eocene of Zagros Basin, Iran[END_REF]

Material and methods

Babaheidar (Babaheydar) is a small town located in Chahar Mahal and Bakhtiari province in west central Iran. It is ca. 50 km from Shahrekord, the capital city of the province (Fig. 1). The area is mostly located in the high Zagros mountain ranges. Numerous marine deposits ranging from late Paleozoic to the end of Cenozoic were found in the area. The 'Konservat Lagerstätte' of Babaheidar was discovered about three decades ago during the construction of a road. The fauna was preliminary reported [START_REF] Jafarian | Preliminary investigation and systematic of Eocene fish from Bakhtiari region[END_REF]), but has not been studied systematically until later [START_REF] Jafarian | Investigation and introduction of some Eocene fish fossils in Zagros Basin[END_REF]. Recent field researches by teams from the University of Isfahan have given thousands of fishes but also plants, crustaceans, insects, and birds [START_REF] Tyler | New genus and species of basal tetraodontoid puffer fish from the Oligocene of Iran, related to the Zignoichthyidae (Tetraodontiformes)[END_REF][START_REF] Garassino | Report on decapod crustaceans from the Eocene of Zagros Basin, Iran[END_REF]. The fossil beds of Babaheidar area belong to the Pabdeh Formation, deposited in the north-western and southern parts of the Zagros Basin mainly during the Eocene and in some areas until the Oligocene and even the Miocene. It mainly consists of white marls, dark to light grey shales rich in planktonic microfauna, and thin to thick-bedded marl limestones (platy limestone) which preserved the macrofauna [START_REF] Garassino | Report on decapod crustaceans from the Eocene of Zagros Basin, Iran[END_REF]. It is about 800 meters thick in the studied area. The Pabdeh Formation is covered by the famous Oligocene/Miocene beds of the Asmari Formation, and it rests on the Eocene terrigenous red deposits of the Kashkan Formation. (Fig. 2).

The sedimentary environment of the Pabdeh Formation based on several studies (microfacies and trace fossil analyses, sedimentary petrology, and micropaleontology) indicates a deep sea environment (upper bathyal zone). The abundance of micrites, glauconite, planktonic foraminifera, fish fauna, and the presence of a pelagic bird support a palaeoenvironment of pelagic marine origin (Paratethys) [START_REF] Aqrawi | The origin of dolomite in the Asmari Formation (Oligocene-Lower Miocene), Dezful Embayment, SW Iran[END_REF][START_REF] Amirshahkarami | Sedimentary facies and sequence stratigraphy of the Asmari Formation at Chaman-Bolbol, Zagros basin, Iran[END_REF][START_REF] Mohseni | Depositional environments and trace fossil assemblages in the Pabdeh Formation (Paleogene), Zagros Basin, Iran[END_REF][START_REF] Garassino | Report on decapod crustaceans from the Eocene of Zagros Basin, Iran[END_REF].

Deposits of the Pabdeh Formation are generally assigned to the Eocene age. However, a study of the nanoplanktons and planktonic foraminifera from one section in the Babaheidar area, in which the pelagic bird was found [START_REF] Peters | Frigidafons babaheydariensis n. sp., ein Sturmvogel aus dem Oligozän des Irans (Aves: Procellariidae)[END_REF], indicated that the studied section is Rupelian (early Oligocene) in age (Hamrsmid & Rogl 2000). Nevertheless, a recent study of the planktonic microfauna from the other sections indicates a middle to late Eocene age [START_REF] Garassino | Report on decapod crustaceans from the Eocene of Zagros Basin, Iran[END_REF]. It is quite possible that the very thick Pabdeh Formation covers a period ranging from the late Eocene to the early Oligocene. The fossil fishes as well as the other fauna are present at several intervals, though they are not abundant in the lower parts.

The fossil insects from Babaheidar area come from several localities. The exact stratigraphic position for some of them is unclear. However, the Orthoptera comes from the upper parts of the succession, early Oligocene in age. It is likely that other insects also come from these intervals because it is the most investigated part of the series for the fossil fishes that are frequently found in these layers. Occurrence. Rupelian, Oligocene, Pabdeh Formation (zone NP 24), Sepid daneh-Kavanak, high Zagros Mountains near the town of Babaheidar (Baba Heydar), western Iran.

Description

Tegmen hyaline, with ca. 11 irregular brown transverse stripes, large, 67.0 mm long, 11.3 mm wide, falcate with apex cut; subcostal area 2.0 mm wide; no visible intercalary longitudinal vein (i.l.v.) in area between R+M+CuA and CuP; RP with ten branches and intercalary longitudinal veins between them; anal veins very long, closely parallel, ending on posterior wing margin ca. 49.0 mm from wing base and 17.3 mm from wing apex; anal area with several rows of very small cells; CuP ending on first anal vein; hindwing large but covered by a tegmen and very poorly preserved, hyaline.

Discussion

The presence of very long anal veins on the tegmen, distally fused with CuPb near posterior wing margin is a character only present in some caeliferan groups (Pauliniidae, Charilaidae, some Catantopidae, and mainly Acrididae) [START_REF] Dirsh | Classification of the acridomorph insects[END_REF]). Among these groups, only some Acrididae have so numerous branches of RP (10 branches) and a similar shape of tegmen elongate, falcate, and with a cut apex (e.g. Cyrtacanthacridinae, Calliptaminae, and Oedipodinae) [START_REF] Dirsh | A preliminary revision of the families and subfamilies of Acridoidea (Orthoptera)[END_REF][START_REF] Dirsh | The African genera of Acridoidea (Orthoptera)[END_REF]. The absence of the intercalary longitudinal vein in the area between R+M+CuA and CuP would exclude affinities with the Oedipodinae. The particular ornamentation of the tegmen with more than 11 transverse dark zones can be found in these subfamilies but also in some Acridinae. The pattern of venation and of ornamentation is very similar to that of the modern Schistocerca gregaria [START_REF] Forskål | Descriptiones Animalium Avium, Amphibiorum, Piscium, Insectorum, Vermium; quae in Itinere Orientall observati Petrus Forskal[END_REF] or Locusta migratoria Linnaeus, 1758, but also of species in the genera Aiolopus [START_REF] Fieber | Synopsis der europäischen Orthoptera. Lotos[END_REF], or Oedaleus Fieber, 1853, etc. (Fard 1990;[START_REF] Garai | Contribution to the knowledge of the Iranian orthopteroid insects I[END_REF]. It is impossible on the sole basis of the wing venation to attribute this fossil to a precise tribe or genus, even if an attribution to the Cyrtacanthacridinae could be the most probable for the great similarity of tegmen with Schistocerca gregaria.

The great similarity of this fossil in size of wings, tegmen shape, venation, and ornamentation with those of a Schistocerca or a Locusta suggests that this insect had flight capacities similar to those of these migratory locusts. Even its long tegmina suggest that it was a long-winged morph able to migrate, unlike the short-winged morph of Locusta migratoria [START_REF] Nishide | Wing dimorphism in the migratory locust, Locusta migratoria: differentiation of wing morph and phase polyphenism[END_REF]. Brown versus green color can have also an impact on the migratory behavior in these locusts [START_REF] Cease | Are color or high rearing density related to migratory polyphenism in the band-winged grasshopper, Oedaleus asiaticus?[END_REF]), but, even if the tegmina of this fossil are colored in brown, it is not possible to be sure that it is the original color of the insect because chemical green colors of insects degrade into brown color during taphonomic process [START_REF] Mcnamara | The taphonomy of colour in fossil insects and feathers[END_REF].

This locust was certainly an active flyer because of the large dimensions and shape of the tegmina. Locusts are known to occur at high altitude during their migration [START_REF] Berland | Données récentes sur le transport aérien d'animaux et de plantes, d'après des recherches faites en avion[END_REF], and to cross over seas and oceans. After its discovery in deep water sediments of the Oligocene Thetys together with its flight capacities, we can infer that it died in open sea and went in deep water sometime after its death because it was obviously dislocated as it only shows fragments of thorax with two forewings and one hindwing. Experiments on modern Orthoptera showed that these insects, when they arrive at the water surface, float for a long time and generally dislocate before sinking [START_REF] Popov | The desert locust (Schistocerca gregaria Forskål) in the Island of Socotra[END_REF]Wilson 1980;[START_REF] Nel | Analyse d'entomofaunes cénozoiques[END_REF][START_REF] Martínez-Delclòs | Insect taphonomy experiments. Their application to the Cretaceous outcrops of lithographic limestones from Spain. Kaupia[END_REF]. The eutrophication in the photic zone of the ancient sea with bloom of nannoplankton [START_REF] Hamrsmid | Biostratigraphy of the Baba Heydar section, Iran[END_REF] favored its sinking and fossilization [START_REF] Nel | Analyse d'entomofaunes cénozoiques[END_REF][START_REF] Paicheler | The aquatic genus Notonecta (Insecta: Heteroptera) as a palaeo-ecological indicator of rhythmite miming sequences in shallow freshwater deposits[END_REF]).

The other insects found in the Pabdeh Formation are small ants (alate males and females Formicidae, 3-5 mm long) (Fig. 5). They were found in the Qanbar Siny section of Babaheidar area. Such ants are flying during their mating swarm. They were lost in the sea. They have also lost body parts (wings, legs, heads or abdomens) and probably stayed a long time on the water surface before sinking.

The concerned area was occupied during the Oligocene by the rather broad Paratethys between the Arabian Peninsula and Central Iran [START_REF] Rögl | Palaeogeographic considerations for Mediterranean and Paratethys seaways (Oligocene to Miocene)[END_REF]

Conclusions

No other record of Iranian fossil insects is known. Further field researches will be necessary to increase our knowledge on the past entomofauna from this country. Such discoveries could be of great interest as direct evidences on the biogeography of this area situated between Africa and Asia. Also, it appears that fossil insects could be more frequent than expected in outcrops of deep marine origin. They could be of great interest to trace ancient behaviors such as intercontinental migrations. UIDGA3-4 (imprint and counter-print); 4. UIDGA5 (imprint); 5. UIDGA7 (imprint); 6.

UIDGA6 (counter-print of UIDGA5); 7-8. UIDGA9-10 (imprint and counter-print); 9.

UIDGA8 (counter-print of UIDGA7); 10-11. UIDGA2 (counter-print of UIDGA1); 12. UIDGA11; 13. UIDGA12; 14. UIDGA13; 15-16 (UIDGA14-15) (imprint and counter-print).

  , none of them was studied. Otherwise only traces of insect activities are putatively known from Iran, viz. Vaez-Javadi (2006: pl. 2, fig. 1) figured fossil leaves with possible arthropod margin traces of activity from the Rhaetian of NE Iran; and Ameri et al. (2014: figs 3-4) figured Middle Jurassic leaves with possible traces of activity on their margins. Dating and palaeoenvironmental conditions for the oldest migrations of locusts between Africa, Europe and Asia are discussed.

  The studied specimens are housed in the Palaeontological Collection of the Department of Geology, University of Isfahan, I.R. Iran (UIDG). The fossils were photographed with a Nikon D40X. We follow the wing venation terminology ofBéthoux and Nel (2002). R, radial vein; RP posterior branch of R; M median vein; CuA anterior branch of cubital vein; CuP posterior branch of cubital vein. Specimen IUDG 16 (two tegmina and a hindwing attached to fragments of the thorax).

  . It is possible to suppose that this locust participated to a migration event between these two lands, similarly to what happens nowadays, covering longer distance than the width of the Oligocene Paratethys between Arabia and Iran (Cressman 1998). The Eocene-Oligocene climate is considered rather dry in the concerned area and more generally in Asia (Dupont-Nivet et al. 2007; Soleimani 2009), fitting well to the presence of migratory locusts living in open grasslands. The cooling and drying of Asian climate at the Eocene-Oligocene boundary caused a turnover of the mammal fauna and flora from a 'perissodactyl-dominant fauna in forests under a warm-temperate climate to small rodent/lagomorph-dominant fauna in a forest-steppe' (Sun et al. 2014). Arabian-Asiatic migrations of locusts could have been initiated and/or favored by these climatic changes during the late Paleogene. More generally it seems that the Acrididae began to diversify only during the Oligocene, maybe also in relation to the development of more open biotas. Fossil Acrididae are very scarce and only known since the late Eocene (Gorochov & Rasnitsyn 2002). The oldest alleged Acrididae, Menatacridium Piton, 1936 from the Paleocene of France, has been discussed by Zeuner (1944) who concluded that its position is unclear. Compression fossils of Acrididae become slightly more frequent in the lacustrine outcrops of the late Oligocene and the Miocene.



  Wilson, M.V.H. 1980. Eocene lake environments: depth and distance-from-shore variation in fish, insect and plant assemblages. Palaeogeography, Palaeoclimatology, Palaeoecology 32:21-44.  Zeuner F.E. 1944. The fossil Acrididae (Orth. Salt.). -Part 4. Acrididae incertae sedis and addendum to Catantopinae. Annals and Magazine of Natural History (11) 11:359-383.

Figure captions Figure 1 .

 captions1 Figure captions

Figure 2 .

 2 Figure 2. Stratigraphic column of the fossiliferous deposits of the Babaheidar lagersttate.

Figure 3 .

 3 Figure 3. Orthoptera Acrididae, specimen IUDG 16. Photograph of imprint. Scale bar: 10 mm.

Figure 4 .

 4 Figure 4. Orthoptera Acrididae, specimen IUDG 16. Photograph of counter-imprint. Scale bar:

Figure 5 .

 5 Figure 5. Hymenoptera: Formicidae (possibly Myrmicinae). 1. UIDGA1 (imprint); 2-3.

Acknowledgements

We sincerely thank Pr Bo Wang and an anonymous referee for their very useful comments on the first version of the paper. MMA thanks Erfan Khosravi (University of Tehran) and Borzu Asgari (Islamic Azad University, Masjed Soleiman) for their support during the field work.

Disclosure statement

No potential conflict of interest was reported by the authors.


Nicolas Belloir

Jérémy Buisson

Olivier Bartheye

Metamodeling NATO Operation Orders: a proof-of

Keywords: Military SoS, Battlefield Engineering, Model-Based Engineering, Operation Orders

published or not. The documents may come

I. INTRODUCTION

Traditionally, military operations are considered as Systemsof-Systems (SoS). Indeed, they usually take place using a combination of complex and social systems: helicopters drop assault sections supported by armored vehicles and artillery for instance. Thanks to interconnected information systems, soldiers can follow and manage on-going operations; in particular, one can mention digital capabilities deployed to the soldier equipped with tablet, communications and digital tools such as connected viewfinders as within the French F ÉLIN system 1 . This evolution was theorized via the concept of networkcentric warfare [START_REF] Alberts | Network centric warfare: developing and leveraging information superiority[END_REF].

These new technological features offer new opportunities to increase the effectiveness of military operations, but require rethinking operational doctrines to effectively integrate them. Mastering these capabilities will impact the way of conducting operations. This will require developing a specific engineering sometimes called battlefield engineering. The latter has to provide tools assisting officers, especially to allow them to take into account the growing complexity of the digital battlefield. One of the challenges will be to adapt the uses so as to rely on information systems to improve efficiency. For example, orders that govern operations will have to be supported by digital technologies, hence leading to a shift from a textual approach in which the order is written in natural language, so difficult to understand automatically, to a format that will make it more easily machine-processable.

On the other hand, software engineering has developed effective methods and tools to take into account the increasing complexity of software development. Some of these tools and methods have been successfully adapted in other areas, such as the Model-Based System Engineering (MBSE) [START_REF] Estefan | Survey of Model-Based Systems Engineering (MBSE) methodologies[END_REF] in system engineering. Indeed, this approach creates a dynamic aiming to shift from textual document-centred design tools and method to model-based ones using graphical languages such as SysML [START_REF]OMG Systems Modeling Language[END_REF] or integrated model-based methods such as Arcadia [START_REF] Roques | Mbse with the arcadia method and the capella tool[END_REF]. The main advantage is clearly to manipulate a single complete model ensuring the inner consistency of the military operation rather than a set of textual documents more or less complete.

In this paper we aim at presenting a proof of concept showing how to manipulate operation orders with a graphical language called OPORD-ML. The latter was created from a metamodel inspired by the NATO standard STANAG 2014 [START_REF]STANAG 2014: Formats for Orders and Designation of Timings, Locations and Boundaries[END_REF] and is tooled by a graphical editor automatically generated from the metamodel using Sirius2 . This makes it possible to write NATO-standard operation orders that can be exported in XMI3 format [START_REF]XML Metadata Interchange (XMI) Specification[END_REF]. The latter is a standard for exchanging metadata information that will easily allow transmitting the written operation order between the various related information systems, or even to autonomous systems such as UAV's for instance.

The remainder of the paper is organized as follows: in the next section we present the military background of this work.

In Section III, we discuss about how Software Engineering approaches could be applied to improve support of new IT capabilities in military operations descriptions. In Section IV, we illustrate the proposed approach through a metamodel, a tool as a proof of concept. Section V discusses some related work while Section VI presents concluding remarks and directions for future work.

II. BACKGROUND

Since the middle of the 1990s, the US Department of Defense has pioneered and elaborated the concept of networkcentric warfare [START_REF] Alberts | Network centric warfare: developing and leveraging information superiority[END_REF]. Acknowledging the rise of information and global digitalization, network-centric warfare denotes attempts to leverage information in the context of military missions and operations. Linking and networking military forces in order to improve shared battlefield awareness provides information dominance and gives advantage in terms of tempo, risk reduction, responsiveness, effectiveness.

While not specifically being a technological issue, networkcentric warfare has driven technology evolution in the military field. Noticeably, command and control (C2) and related functions can be supported by dedicated information systems. For instance, US Army's Command Post of the Future (CPOF) provides commanders with a tool to manage views of the battlefield for the mission they command, that is, produce, visualize and share any battle-related artifacts. By means of a communications network, CPOF bidirectionally exchanges data, e.g., with the Tactical Ground Reporting (TIGR) system equipping warfighters at the company level and below, down to soldiers. Likewise, the French Army has acquired similar technologies, for similar purpose, such as Force Command Information System (SICF), Regiment Information System (SIR) and Terminal Information System (SIT), from the command post down to the section level. The Integrated Infantryman Equipment and Communications (F ÉLIN) equipment provides soldiers with a portable electronic platform that can be used in conjunction with SIT by means of standard PR4G radio, and the Félin Information Network (RIF) allows information sharing within the infantry section. The more recent Scorpion Combat Information System (SICS) ensures exchange of information between different levels of command and across arms.

Such battle management systems currently provide several features. Blue force tracking consists in sharing the tactical situation amongst warfighters. The command post and warfighters follow the progress of friendly forces during operation execution. Red force tracking aggregates situation reports and intelligence in order to achieve the same with respect to enemy forces. Besides coordination of friendly forces, having such situation awareness helps reducing the fog of war uncertainty when making tactical decision during operation execution, in order to better identify and exploit opportunities, and gain tactical initiative.

Having a consistent information system also allows connecting target acquisition systems to weapons systems. For instance, reconnaissance units may interact with artillery units in the back of the battlefield in order to designate targets. Such a connection may even allow front units to trigger fire of back units. Front units may indeed detect and strike threats in advance, before weapon units located more in the back can detect them.

Current battle management systems also offer tools and interfaces for the elaboration of orders, as well as the network for their transmission. Depending on the tradition of each nation's or coalition's army, a framework depicts how orders are structured, organized, formatted. In the following, we use NATO's STANAG 2014 [START_REF]STANAG 2014: Formats for Orders and Designation of Timings, Locations and Boundaries[END_REF] as an illustration, because it is an unclassified document. STANAG 2014 makes the distinction between three kinds of orders. First, warning orders (WngO) tell units essential information to initiate their preparations in prevision of a forthcoming mission. Then, operation plans (OPLAN) and operation orders (OPORD) describe the operation itself. At the plan phase, the OPLAN aims at preparing a future operation based on an anticipated situation. When the conditions are met and operation execution is determined, the OPLAN turns into an OPORD. Fragmentary orders (FRAGO) are issued either when a partial content of the order is only needed or when a preexisting OPORD has to be amended quickly, without spending time in the elaboration of a complete new OPORD.

Regardless of its kind, any STANAG 2014 order is structured into 5 well-defined sections. A situation section first introduces the order mainly with a description of enemy and friendly forces. Then a mission section states the objectives of the operation. An execution section explains the intent, concept and manoeuvre of the operation in order to fulfill the mission. An administration / logistics section describes operation support activities. Last, a Command and signal section describes the communication and electronic policy of the operation. Additional concerns, such as intelligence, fire support, engineer, may be addressed by domain-specific annexes attached to the order. In addition to the textual description of the order, an overlay order (OvO) may provide a graphical representation of forces and manoeuvre using standard symbols from NATO's APP-06 [START_REF]APP-06 NATO Joint Military Symbology[END_REF] drawn a on map. WngO, OPLAN, OPORD, FRAGO and OvO are related to one another when describing a single operation.

In this context, orders are artifacts of communication along the chain of command. It is used in the overall workflow of operation planning and execution, and battle management systems support this workflow. After sending his order to his subordinate, the higher commander briefs this subordinate to explain his order and its underlying philosophy. Then, during the back-brief phase, the subordinate explains to his commander his understanding of the operation, in order to ensure it is correct. Analysis and reasoning lead the subordinate to elaborate his own orders in order to fulfill the objective assigned by the commander, yielding to lower-level orders. Then the order is played during mission rehearsal to ensure that the desired effect is actually obtained without unanticipated drawbacks, before operation execution. For each order at a given hierarchical level in the chain of command, several orders are derived at the level below, down to the lowest level, that implement the operation prescribed by the order.

However, a battle management system is not only about having an information system that stores and broadcasts orders along the chain of command. Like described by Khimeche [START_REF] Khimeche | C2-Simulation interoperability benefits[END_REF], a battle management system supporting C2 may be connected with simulators that provide added value. Simulation during briefing and back-brief help improve mutual understanding by providing visual animated illustrations. Simulated rehearsal is a testing step that helps detecting potential weaknesses and issues in the forthcoming manoeuvre. During analysis and reasoning, simulation of the course of action augment the warfighter's expertise with computational analysis and comparison for arbitration. During execution, simulation helps detect how much the situation diverges from the planned operation, in order to help decide whether issuing a FRAGO is relevant. More than an information system, the battle management system is a backbone on which tools are pluggedin in order to assist the operation.

Analysis and design of a manoeuvre itself is subject to methodology. Yakovleff's book [START_REF] Yakovleff | Tactique théorique[END_REF] is one of the contemporary references in the French army. Beyond its theoretical aspect, Yakovleff's framework is, of course, related to the content of an order, hence structuring and tactically consistent design activities in terms of space, enemy, time. The manoeuvre is built of elementary actions. Based on historical experience, a collection of well-identified tactical moves hints typical patterns of actions depending on the context. At a more generic level, doctrines prescribe the overall approaches in terms of tactical organization and actions in response to anticipated risks, based on acceptance by political authorities.

In summary, elaboration of military operation orders is similar to an engineering process. The underlying workflow encompasses expression of mission by the commander, its transmission to subordinate, elaboration of the manoeuvre, and its validation before execution. The battle management system provides a collaboration bus that support the workflow. On the one hand, its provides a communication system between stakeholders during the planning phase, for the elaboration of orders. Simulators are connected to this bus in order to assist in this task. On the other hand, blue and red force tracking provides supervision tools that help conducting during operation execution.

III. VISION

We envision that battle management systems like described in previous Section II are still in early stage. Beyond current communication systems targeted at the exchange of information along the command chain and between the battlefield and commander, battlefield digitalization enables high-level tools that can improve the management of a military operation, when seen as a System-of-Systems. Such high-level tools are desirable in order to improve efficiency of the management of the operation. Namely, information shall be propagated and verified in order to ensure overall consistency.

Relying on a model, that is, complete information based on a structured format, rather than informal text or graphical notations based on natural language, allows to build on practices in other fields such as system engineering, requirement engineering and software engineering. We indeed observe similarities between these fields and military operation engineering: the mission describing operation objectives is similar to elicited requirements; doctrines describing typical mission organizations are similar to architectural styles and patterns; tactical moves are similar to design patterns; the overall military operation involves units in order to fulfill the assigned mission like a System-of-Systems involves subsystems to achieve requirements.

In the fields of system engineering, requirement engineering and software engineering, models have proved to be a valuable artifact to describe system architecture, to elicit requirements and to document system implementation while maintaining overall traceability. Holistic traceability, relating doctrines, tactical moves, operations at all levels in the chain of commands, live tracking during operation execution, and lessons, all in a single model, will provide comprehensive description of an operation and related information. Model-Driven Engineering (MDE) enables generation of consistent infrastructure for tools from the metamodel describing the abstract language of models. In our case, in addition to a graphical editor sketched in Section IV-C, for instance, such tools may include automatic generation of textual orders and of overlay orders from the model, and automatic propagation of information along the chain of command, as well as artifacts ensuring interoperability with simulators. We expect that configuration and instructions for battlefield UAVs and upcoming robots will derive from the model as well. Relating order elements to higher-level orders and to doctrines not only enables documentation and explainability of commander's decision, it will also help setting up analysis and verification in order to ensure overall consistency with the operation framework.

Our vision is that model-based technologies and practices in system, requirement and software engineering pave the way towards achieving the promises of network-centric warfare.

IV. METAMODELIZATION OF NATO ORDERS A. Model-Based System Engineering

Model-Based System Engineering (MBSE) is a system engineering practice using the capabilities of computer technologies to describe, through models, concepts, and languages, both the problem and its solution. It is recognized as a promising direction for dealing with system and SoS complexity [START_REF]Systems Engineering Vision 2020, version 2.03[END_REF]. One of the MBSE core ideas is to store into models the engineering knowledge of a specific domain. Domain semantics is stored in a metamodel. Abstract and concrete syntaxes are defined from the metamodel insuring conformity to the described semantics. Transformations are another MDSE core principle: transforming a model conforming to a metamodel M1 into another model conforming to metamodel M2 allows

B. A metamodel for NATO operation orders

We present here a metamodel inspired by NATO's STANAG 2014 [START_REF]STANAG 2014: Formats for Orders and Designation of Timings, Locations and Boundaries[END_REF] norm. The metamodel is still under development and is not yet fully compliant with the norm. Thus, it must be considered together with the generated tool (presented in the next subsection) as a proof of concept.

Figure 1 shows the proposed metamodel. The Order metaclass is the core of the metamodel. An Order is constituted by the following metaclasses: Situation, Mission, Execution, Official, Distribution, Security Classification, Acknowledge, Com-mand and Signal, Administration and Logistic and Annexes. The Situation metaclass describes the context of the operation. The latter is compound by Weather and Terrain metaclasses. Moreover it helps to describe forces to be considered: three sorts of Units can be expressed using the Enemy, Allied and Attachments and Detachments metaclasses. Units are characterized by their name and by a textual description (name and information meta-attribute). If more precision needed, it is possible to use the inheritance mechanism to specialise each unit. For instance an unit can be an Infantry Unit. We illustrate this possibility in our metamodel by proposing 4 specializations: Armored Unit, Mechanized Unit, Custom Unit and Infantry Unit. It will be extended to encompass all the possibili-ties. The Execution metaclass is used to describe sequences of actions that must be realized in the context of the designed operation. We limited possible actions to combat and manoeuvre identified by Tasks and Missions to Combat Support Units and Tasks and Missions to Combat Manoeuvre Units metaclasses. Both can be coordination actions using Coordinating instructions. Both are related to one or more friendly or enemy units by links to Unit metaclass. OCL rules prevent being able to describe a combat action between two friendly units in order to avoid friendly fires, to illustrate consistency rules that can be encoded in the metamodel. We define a domain-specific language called OPORD-ML and based on the previous metamodel. It allows to describe operation orders. A tool implementing OPORD-ML was generated using Sirius. It defines a graphical concrete syntax implementing the proposed metamodel. This tool is generated using MDE principles. The user, here the commander, builds the order he wants to design using the Palette illustrated by Figure 3. Figure 2 shows part of an order designed using OPORD-ML and using the generated tool. It is decomposed in two parts. The first one concerns the original situation. All the mandatory data are available: weather, terrain, attachments and detachments, enemy and allied forces. The second part concerns the description of the operation in terms of execution phases. In this example, three sequential phases were described. Arrows specify the dependencies. For instance, Phase 1 -Attack concerns one allied unit and one enemy unit. Semantics embedded in the metamodel controls the order: it makes impossible to attach an allied unit as a target for an other allied unit in a combat action.

C. OPORD-ML modeling language and tool

Each model element is characterized by a property window. For instance, Figure 4 shows that weather is given by Min Fig. 6. Attack Action Window temp, Max temp, Forecast description and Date. Figure 5 shows the properties of an allied unit. Properties windows depend on the metamodel. It can be more complex as shown by Figure 6. Friendly and Enemy units are automatically added or deleted using ad-hoc buttons. Thus it is not possible to add a unit that was not identified in the situation part of the order.

V. RELATED WORKS

Early work in the area [START_REF] Salisbury | Command And Control Simulation Interface Language (CCSIL): Status Update[END_REF], [START_REF] Ogren | EAGLE Combat Model Battle Management Language (BML)[END_REF] originally aimed at interoperability with simulators, identified the need for a shift from free text towards highly structured data.

The Battle Management Language (BML) [START_REF] Carey | Standardizing battle management language -a vital move towards the army transformation[END_REF], [START_REF]Standardizing Battle Management Language -Facilitating Coalition Interoperability[END_REF] is an early attempt to define a language to represent orders. In this work, Carey et al. acknowledge the need for an interchange format that allows command-and-control systems, simulators and upcoming robotic forces to exchange orders and share commun operational picture. Instead of relying on NATO's standard format like we do, Carey et al. based their format on the 5Ws framework: Who, What, When, Why, Where. As a result, their language is centered on the task concept, and the operation, for instance, is not reified per se. To some extent, BML and C-BML allow only the description of tasks within the Execution object of our language. Regarding the methodology, Carey et al. define their language in terms of a logical model, with relational data base in mind. Our approach based on models offers a richer description of the language, and technologies like, e.g., Neo4EMF [START_REF] Benelallam | Neo4emf, A Scalable Persistence Layer for EMF Models[END_REF] have shown that model-based engineering is flexible enough to use various storage backends, including relational and graph-based data bases in addition to the XMI format.

The Coalition Battle Management Language (C-BML) proposed in [START_REF] Sudnikovich | Extensible battle management language as a transformation enabler[END_REF] and accepted in [START_REF] Blais | Coalition-battle management language (c-bml) study group final report[END_REF] extends the work of Carey et al. with service oriented architecture, web service and XML technologies. Gustavsson et al. [START_REF] Gustavsson | Object-Orientated implementation of Grammar Based Battle Management Languages[END_REF] propose to derive an object-oriented model from the C-BML grammar description. They do not improve over BML with respect to our observations, especially with respect to conformance with NATO STANAG 2014.

VI. CONCLUSION

To gain the benefits of augmenting the digital aspects of SoS, in this paper we argue that new tools and methods need to be developed to allow for their integration into advanced doctrines. We proposed the OPORD-ML modeling language allowing to graphically write an operation order. It is based on a metamodel ensuring tactical consistency and supported by an automatically generated tool. This language is inspired from the standard NATO STANAG 2014 [START_REF]STANAG 2014: Formats for Orders and Designation of Timings, Locations and Boundaries[END_REF]. Thus, operation orders can be exported in the XMI standard format [START_REF]XML Metadata Interchange (XMI) Specification[END_REF], allowing to be transferred and interpreted automatically.

Currently, our work must be seen as a proof-of-concept. OPORD-ML is fully compliant with the standard NATO STANAG 2014. We are working on its extension in order to correct that. Moreover, operation orders are generally created at a high hierarchical level. They are then refined iteratively at each lower hierarchical level. We plan to identify decomposition mechanisms to realize automatically part of the refinement. Another topic of interest would be to extend the metamodel in order to link it with geographical concepts. Indeed, operation orders are always redacted considering a map of the targeted battlefield, which has to be integrated within the model.

Fig. 1 .

 1 Fig. 1. Nato OPORD inspired Metamodel

Fig. 2 .

 2 Fig. 2. OPORD

Fig. 4 .

 4 Fig. 4. Weather Properties Window

Integrated Infantryman Equipment and Communications: https:// www.safran-group.com/defense/soldier-modernization

Sirius: http://www.eclipse.org/sirius/

XMI: XML Metadata Interchange

Acknowledge Official