Diego Carnevale
email: diego.carnevale@ens.fr

Philippe Pelupessy

Geoffrey Bodenhausen

Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel?

Keywords: Orientational Inequivalence, Quantum Tunneling, Magic-Angle Spinning, Solid-State NMR, Anisotropic Motional Dynamics, Chemical Exchange, Homonuclear Dipolar Interaction, Chemical Shift Anisotropy, Cross-Term Splitting

Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel? 2

Introduction

Despite its deceptive simplicity, H2O shows a variety of surprising properties, whether in bulk liquids or solids, confined in nanocavities, or near the surface of macromolecules such as proteins or nucleic acids. [START_REF] Ball | Water -an enduring mystery[END_REF][START_REF] Finney | Water? What's so special about it? Phil[END_REF][START_REF] Angell | Insights into Phases of liquid water from study of its unusual glass-forming properties[END_REF][START_REF] Israelachvili | Role of hydration and water structure in biological and colloidal interactions[END_REF][START_REF] Algara-Siller | Square ice in graphene nanocapillaries[END_REF][START_REF] Koga | Formation of ordered ice nanotubes inside carbon nanotubes[END_REF][START_REF] Kolesnikov | Anomalously soft dynamics of water in a nanotube: a revelation of nanoscale confinement[END_REF] To date, there is no established general consensus on the structural organization of H2O molecules on the atomic scale. Hydrogen bonds and dynamic processes such as proton exchange further complicate the scenarios. Quantum effects have been considered to shed light on this molecular jumble. [START_REF] Hernández De La Peña | Temperature dependence of quantum effects in liquid Water[END_REF][START_REF] Paesani | Quantum effects in liquid water from an ab initio-based polarizable force field[END_REF][START_REF] Drechsel-Grau | Collective proton transfer in ordinary ice: local environments, temperature dependence and deuteration effects[END_REF] Recently, quantum tunneling dynamics have been reported for water molecules confined in nanocavities of beryl crystals [START_REF] Gorshunov | Quantum behavior of water molecules confined to nanocavities in gemstones[END_REF][START_REF] Kolesnikov | Quantum tunneling of water in beryl: a new state of the water molecule[END_REF] and in gasphase clusters like (H2O)6. [START_REF] Richardson | Concerted hydrogen-bond breaking by quantum tunneling in the water hexamer prism[END_REF] The chemical shift anisotropy (CSA) of protons in water molecules is of particular interest since it has the potential of lifting the magnetic equivalence of the two protons. Inter alia, this may give access to para-water, i.e., by inducing an excess or a deficiency of the population of the singlet state    =+  , either by going to very low spin temperatures, [START_REF] Mammoli | Challenges in preparing, preserving and detecting para-water in bulk: overcoming proton exchange and other hurdles[END_REF][START_REF] Meier | Spin-isomer conversion of water at room temperature and quantum-rotor-induced nuclear polarization in the water-endofullerene H2O@C60[END_REF] or by suitable manipulations with radio-frequency (rf) pulses. [START_REF] Eisendrath | NMR of protons in gypsum. I. Experimental proof of the existence of four thermodynamic invariants[END_REF][START_REF] Eisendrath | NMR of protons in gypsum. II. Independent observation of molecular reorientations and molecular diffusion[END_REF] In the absence of motional averaging, the shielding tensors of the two protons in a water molecule have different orientations with respect to one another. Except if the crystal environment breaks the symmetry of the H2O molecules (which is not the case in our study), the two protons are chemically equivalent, i.e., have identical tensors in their principal axis systems, with identical isotropic shifts iso, identical anisotropies CS and identical asymmetry parameters ηCS. When immersed in a magnetic field however, the difference of the orientations of the two shielding tensors renders the two protons inequivalent, except for some specific orientations.

More specifically, although with identical isotropic shifts iso and identical rank-2 anisotropic shift tensors  

    =   =
, where the principal axis frame P of spin I is

Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel?

4
arbitrarily chosen as common frame C. We refer to this type of non-equivalence as orientational inequivalence.

We have investigated the protons of water molecules trapped in barium chlorate monohydrate, Ba(ClO3)2•H2O, by means of variable temperature magic-angle spinning (MAS)

proton NMR spectroscopy over a temperature range 110  T  300 K. This sample bears analogies to hydrated gypsum CaSO4•2H2O that features characteristic proton NMR spectra that were first described by Pake. [START_REF] Pake | Nuclear resonance absorption in hydrated crystals: fine structure of the proton line[END_REF] Similar dynamic processes have been studied in Ba(ClO3)2•D2O by

Long et al. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] by means of deuterium NMR spectroscopy, albeit over a narrower temperature range 154  T  291 K. We have independently simulated the deuterium lineshapes observed by Long et al. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] and found good agreement with their exchange rates. Proton and oxygen-17 NMR investigations of the same sample have also been reported. [START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF][START_REF] Keeler | 17 O NMR investigation of water structure and dynamics[END_REF] In this work, we show numerical simulations of proton spectra over an extended temperature range 110  T  300 K that take into account the presence of exchange and agree remarkably well with the unusual sideband patterns of our experimental MAS spectra. The experimental evidence is compared with previous 2 H and [START_REF] Ball | Water -an enduring mystery[END_REF] H NMR studies. Within experimental errors over the 140  T  190 K temperature range, we find an activation energy which is higher than that reported in the literature. This finding is corroborated by density functional theory (DFT) calculations of the thermal activation energy for the two-fold flip about the H-O-H bisector. Deviations from linearity in the Arrhenius plot are rationalized by means of the Bell-Limbach model [START_REF] Bell | The Tunnel Effect in Chemistry[END_REF][START_REF] Buntkowsky | H-solid state NMR studies of tunneling phenomena and isotope effects in transition metal dihydrides[END_REF] that takes into account incoherent quantum tunneling contributions at low temperatures. The unusual spinning sideband patterns observed at low temperatures result from CSA-dipole cross terms that lead to fourth-rank interactions that cannot be eliminated by spinning about the magic angle. [START_REF] Maricq | NMR in rotating solid[END_REF][START_REF] Ashbrook | Second-order cross-term interactions in high-resolution MAS NMR of quadrupolar nuclei[END_REF]

Results and discussion

The crystal structure of Ba(ClO3)2•H2O is shown in involved in two hydrogen bonds of length 1.92 Å with oxygen atoms of two different chlorate anions. The oxygen atom of each water molecule is coordinated with one Ba 2+ cation, at a distance of 2.74 Å, along a vector parallel to the b-axis. This structure defines channels of ca. 4.5 Å across that run parallel to the c-axis, where water molecules are found at intervals of 4.88 Å, with a 180° alternation of the orientations of the H-O-H bisecting vectors. Intermolecular dipolar couplings between the closest protons ( 2 kHz) are therefore much smaller than intramolecular couplings (ca. 30 kHz), so that each water molecule may be considered to be reasonably isolated. This is consistent with the quasi-ideal Pake patterns observed in static proton NMR spectra of this sample. [START_REF] Carnevale | Solid-state NMR measurements and DFT calculations of the magnetic shielding tensors of protons of water trapped in barium chlorate monohydrate[END_REF][START_REF] Carnevale | Exciting wide NMR spectra of static solid samples with weak radiofrequency fields[END_REF] These patterns are compatible with the fact that the dominant intramolecular dipolar interactions is not affected by the swapping of the protons through a two-fold flip about the H-O-H bisector or through tunneling. Therefore, the intramolecular dipolar Hamiltonian is an invariant of motion in our case. [START_REF] Eisendrath | NMR of protons in gypsum. I. Experimental proof of the existence of four thermodynamic invariants[END_REF][START_REF] Eisendrath | NMR of protons in gypsum. II. Independent observation of molecular reorientations and molecular diffusion[END_REF] On the other hand, fast swapping of the protons results in a partial averaging of the two shielding tensors, resulting in two collinear tensors, thus removing their orientational inequivalence above 200 K. In this Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel?

6 case, we can speak of a quasi-invariant of motion for the CSA Hamiltonian. We have previously investigated the CSA tensors of the protons of trapped water molecules by solid-state NMR under both static and magic-angle spinning conditions at 300 K. [START_REF] Carnevale | Solid-state NMR measurements and DFT calculations of the magnetic shielding tensors of protons of water trapped in barium chlorate monohydrate[END_REF] The values measured at room temperature agree remarkably well with those calculated by means of DFT methods for a periodic lattice when fast 180° flips about the H-O-H bisectors are assumed to occur rapidly. In the absence of motional averaging, these calculations show that the two protons have shielding tensors that are identical in size and shape. [START_REF] Carnevale | Solid-state NMR measurements and DFT calculations of the magnetic shielding tensors of protons of water trapped in barium chlorate monohydrate[END_REF] In this study, we have repeated these DFT calculations at a higher level of theory. The resulting NMR parameters are CS = -15.4 ppm, CS  = 0.1, with Euler angles relating the two shielding tensors (90.13°, 65.37°, 90.12°), in good agreement with our previous study. [START_REF] Carnevale | Solid-state NMR measurements and DFT calculations of the magnetic shielding tensors of protons of water trapped in barium chlorate monohydrate[END_REF] The angle  = 65.37° between the two main components zz of the shielding tensors is compatible with the tetrahedral angle, i.e., ' = (180-65.37)° = 114.63°. Figure 2 shows numerical simulations of proton spectra assuming these latter inhomogeneities and some degree of disorder in the physical sample. These simulations show a remarkable agreement with the experimental low temperature lineshapes. The exchange rates kex were chosen to obtain the best resemblance with the experimental spectra and by minimizing differences in linewidths. We ascribe subtle discrepancies between relative intensities of spinning sidebands in the experimental and simulated spectra to the fact that T1 depends on the orientation of the crystallites and is therefore inhomogeneous across the spinning sideband manifolds. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF][START_REF] Vold | Deuterium relaxation in molecular solids[END_REF] Figure 4 shows the experimental and simulated linewidths (black and red data points) extracted from the lineshapes of Fig. 3 vs 1/T. The overall agreement is very satisfying. The coalescence region, characterized by a marked increase in linewidths, is clearly identified. A 900

Hz line broadening allows one to match the FWHH at high temperatures. However, a mismatch Full widths at half height (FWHH) of the experimental MAS spectra of the left panel of Fig. 3 (black data points). Analogous FWHH of the simulations of the right panel of Fig. 3 (red data). Values were calculated with respect to the centerband. For T  140, the average FWHH for the two lines of the 'doublet' is reported. in linewidths is observed in the low-temperature region. We ascribe this discrepancy to errors in the Euler angles calculated by DFT methods, despite the improved accuracy of recent calculations. Indeed, a variation of these Euler angles leads to a fine structure (which is perceived as broadening) of each component of the low-temperature doublets (see Supporting Information, Fig. S1). It is also worth noting that the lack of temperature dependence of the linewidths in the first three and last three experimental and simulated spectra indicate that variations kex rates do not affect the linewidths in these ranges. This is in agreement with the observation that the three spectra in the high-temperature range 200  T  300 K of Fig. 3 do not show any significant changes in lineshape. Similarly, the experimental proton spectra show little variation in the low-temperature range 110  T  130 K. dots. The error bar for the i th data point lnkex (i) was set to lnkex (i+1) -lnkex (i-1) to represent the worst reasonable uncertainties. Larger discrepancies in lnkex (i) would require to alter the order of the simulated spectra of Fig. 3 with respect to the experimental ones. A linear fit (dashed line) yields an activation energy Ea = 36.0  1.5 kJmol -1 and pre-exponential factor A = 8.4  10 15 s -1 . and A = 1.9 × 10 14 ). In grey, these latter data points are adapted for 1 H by assuming the usual kinetic isotope effect (Ea( 1 H) = Ea( 2 H)/2 1/2 = 21.7 kJmol -1 , keeping the same pre-exponential factor), in agreement with Ea( 1 H) = = 19.2 kJmol -1 estimated from the temperature dependence of T1( 1 H) by Silvidi. [START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF] The red line indicates an Arrhenius equation with a slope or thermal activation energy Ea = 44.2 kJmol -1 calculated by planewave pseudopotential DFT methods with an intercept or pre-exponential factor A = 1.4 × 10 18 chosen so as to match the kex of the 1 H case (grey data point) at high temperature (T = 300 K, left-most grey data point). The green line has the same slope as the red one, with a pre-exponential factor equal to that of Long et al. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel?

11

The cyan data points represent the 2 H exchange rates as measured by Long et al. which lead to an activation energy Ea = 30.4 kJmol -1 and pre-exponential factor A = 1.5  10 14 s -1 . [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] In grey, these latter data points are adapted for the 1 H case by assuming the usual kinetic isotope effect (Ea( 1 H) = Ea( 2 H)/2 1/2 = 21.7 kJmol -1 , keeping the same pre-exponential factor), [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] in agreement with Ea( 1 H) = 19.2 kJmol -1 estimated by Silvidi from the temperature dependence of T1( 1 H). [START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF] Therefore, over the range 140-190 K, our measurements result in exchange rates that are much closer to those found for 2 H by Long et al. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] than to those estimated for 1 H by Silvidi. [START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF] Furthermore, within our experimental errors, the activation energy measured in our [START_REF] Ball | Water -an enduring mystery[END_REF] H study is comparable if not larger than that found for 2 H. In order to rationalize these findings, we investigated the activation energy associated with a classical thermally activated exchange process by means of a series of single-point energy (SPE) calculations performed with planewave pseudopotential DFT methods on the periodic lattice. A single water molecule within one geometry-optimized unit cell was rotated about the H-O-H bisector in steps of 18° from 0 to 180°. SPE calculations were performed on the resulting geometries. The corresponding energy profile, shown in Supporting Information (Fig. S2), yields an activation energy Ea = 44.2 kJmol -1 . The red line in Figure 5(a) indicates the slope obtained with this calculated value when the pre-exponential factor A = 1.4  10 18 s -1 is chosen to match the experimental rate kex = 3  10 10 s -1 for 1 H at T = 300 K. In contrast, the green line indicates an Arrhenius plot obtained for Ea = 44.2 kJmol -1 estimated by DFT methods in this study but keeping the same pre-exponential factor A = 1.9 × 10 14 found for 2 H by Long et al. [START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] In this latter case the agreement with the experimental data points is evidently very poor. Both red and green lines in Fig. 5(a) indicate that both activation energies and pre-exponential factors determined experimentally in this study and those reported in the literature are smaller than those calculated by DFT methods.

In order to explain these discrepancies, we considered effects of quantum tunneling, as described by Bell. [START_REF] Bell | The Tunnel Effect in Chemistry[END_REF] More generally, quantum tunneling phenomena have been investigated in a variety of chemical contexts, with particular attention to their role in kinetic isotope effects. [START_REF] Shiner | The Arrhenius parameters of the deuterium isotope rate effect in a basepromoted elimination reaction 1 : evidence for proton tunnelling[END_REF][START_REF] Cha | Hydrogen tunneling in enzyme reactions[END_REF][START_REF] Limbach | Primary kinetic HH/HD/DH/DD isotope effects and proton tunnelling in double proton-transfer reactions[END_REF][START_REF] Limbach | NMR and INS Line Shapes of Transition Metal Hydrides in the Presence of Coherent and Incoherent Dihydrogen Exchange[END_REF][START_REF] Wehrmann | Coherent D2 rotational tunneling and incoherent D2 dynamics in a solid non-classical RuD2 complex studied by 2H solid state NMR spectroscopy[END_REF][START_REF] Limbach | Arrhenius curves of hydrogen transfers: tunnel effects, isotope effects and effects of pre-equilibria[END_REF][START_REF] Meisner | Atom tunneling in chemistry[END_REF][START_REF] Horsewill | Quantum tunnelling aspects of methyl group rotation studied by NMR[END_REF][START_REF] Pohl | Proton transfer in guanine-cytosine base pair analogues studied by NMR spectroscopy and PIMD simulations[END_REF] Within this framework, the probability of a quantum object such as a proton with an energy smaller than the potential barrier to appear on the other side is taken into account in the evaluation of the apparent rate constant of the process. This probability is greater for particles Cross-term splittings due to the orientational inequivalence of proton magnetic shielding tensors: Do water molecules trapped in crystals hop or tunnel?

12
with smaller masses and for narrower barriers. Bell's model has been extensively utilized by

Limbach and collaborators for a variety of systems. [START_REF] Buntkowsky | H-solid state NMR studies of tunneling phenomena and isotope effects in transition metal dihydrides[END_REF][START_REF] Limbach | Primary kinetic HH/HD/DH/DD isotope effects and proton tunnelling in double proton-transfer reactions[END_REF] The blue line in Fig. 5 the orientational inequivalence of the two shielding tensors was removed by assuming that the two shielding tensors were collinear (and orthogonal to the dipolar tensor). This assumption leads to a single spinning sideband manifold, with an overall spread that reflects the strength of the intramolecular dipolar-dipole coupling, like the width of the Pake pattern of a static nonspinning powder. This proves that the cross-term splitting between the two sideband manifolds observed at low temperatures requires orientational inequivalence. Furthermore, Fig. S3(c) shows a simulated spectrum where the dipolar interaction was removed. One observes again a single spinning sideband pattern, albeit with a much narrower range because the CSA at 18.8 T is much smaller than the dipolar interaction. Clearly, the cross-term splitting between the two interleaved sideband manifolds observed in proton MAS spectra at low temperatures only occur when both above conditions are fulfilled simultaneously.

The unusual splitting observed at low temperatures results from cross terms between the predominant intramolecular homonuclear dipolar coupling and the two shielding tensors provided these are orientationally inequivalent. Similar cross terms have been extensively studied in solution-state NMR, where they lead to cross-correlated relaxation effects. [START_REF] Kumar | Cross-correlations in NMR[END_REF] In solids, the resulting cross terms lead to fourth-rank effects that cannot be eliminated by magic angle spinning. [START_REF] Maricq | NMR in rotating solid[END_REF][START_REF] Ashbrook | Second-order cross-term interactions in high-resolution MAS NMR of quadrupolar nuclei[END_REF] This is shown in Figure 6, where the effect of the spinning rate on the cross-term splittings observed in this study is investigated by numerical simulations.

Experimental attempts at either saturating or inverting one of the two spinning sideband manifolds observed at low temperature by means of a train of rotor-synchronized DANTE pulses [START_REF] Caravatti | Selective pulse experiments in high-resolution solid state NMR[END_REF][START_REF] Carnevale | Broadband excitation in solid-state NMR of paramagnetic samples using Delays Alternating with Nutation for Tailored Excitation ('Para-DANTE[END_REF][START_REF] Carnevale | Solid-state proton NMR of paramagnetic metal complexes: DANTE spin echoes for selective excitation in inhomogeneously broadened lines[END_REF][START_REF] Lu | Broadband excitation in solid-state NMR using interleaved DANTE pulse trains with N pulses per rotor period[END_REF] were unsuccessful, since the other non-irradiated manifold was also affected (data not shown). This can be understood with the help of the simulations shown in Supporting Information (Fig. S4 andS5), where either single-spin or single-transition operators were detected in both static and spinning cases.

Conclusions

The analysis of hitherto unreported cross-term splittings between two manifolds of spinning sidebands in proton MAS spectra has allowed us to investigate the dynamics of protons

Experimental and computational details

The sample was purchased from Sigma-Aldrich and kept overnight in an oven at 70°C

prior to use in order to minimize adsorption water. All proton MAS spectra were recorded at the Ecole Normale Supérieure on a wide-bore Bruker 800 spectrometer (18.8 T) with an Avance-III console, using 3.2 mm ZrO2 rotors spinning at 10 kHz in a triple-resonance low-temperature MAS probe designed for DNP experiments, but without microwave irradiation. The rf-field amplitude of the proton pulses was 1 = 125 kHz, corresponding to a pulse length τp = 2 μs for a 90° nutation angle. All spectra were acquired with rotor-synchronized solid echoes (90° - -90° - -detection.) Exorcycle [START_REF] Bodenhausen | Suppression of artifacts in two-dimensional J spectroscopy[END_REF] was applied to the second 90° pulse. DFT calculations were conducted with Quantum ESPRESSO 6.3. [START_REF] Giannozzi | QUANTUM a modular and open-source software project for quantum simulations of materials[END_REF][START_REF] Giannozzi | Advanced capabilities for materials modelling with Quantum ESPRESSO[END_REF] Geometry optimization was performed on the known crystal structure [START_REF] Sikka | Neutron-diffraction refinement of the crystal structure of barium chlorate monohydrate[END_REF] NMR calculations were performed on the crystal structure optimized with a cut-off energy of 100 Ry, 4  4  4 k-point grid and norm-conserving pseudopotentials for all atoms. The gaugeincluding projector augmented wave (GIPAW) algorithm was employed to reconstruct allelectron wave functions in the presence of a magnetic field. [START_REF] Pickard | All-electron magnetic response with pseudopotentials: NMR chemical shifts[END_REF][START_REF] Yates | Calculation of NMR chemical shifts for extended systems using ultrasoft pseudopotentials[END_REF] The graphical tool MagresView [START_REF] Sturniolo | Visualization and processing of computed solid-state NMR parameters: MagresView and MagresPython[END_REF] was utilized to extract the sets of Euler angles relating the orientations of the relevant tensorial interactions when expressed in a common crystal frame by a DFT calculation of NMR parameters.

The numerical simulations in Fig. 3 of MAS spectra obtained by Fourier transformation of solid echoes with inclusion of chemical exchange were carried out by means of the master equation, as applied to solids by Duer and Levitt: [START_REF] Ernst | Principles of Nuclear Magnetic Resonance in One and Two Dimensions[END_REF][START_REF] Duer | Time-domain calculation of chemical exchange effects in the NMR spectra of rotating solids[END_REF] 𝝆(𝑡) = exp{(-𝑖𝐿 𝑐𝑜ℎ + 𝐿 𝑒𝑥 )∆𝑡}𝝆(0),

where Lcoh is the commutation superoperator describing coherent spin dynamics due to the anisotropic chemical shift and dipolar Hamiltonians, i.e., all interactions relevant in our case as far as the water molecules can be considered to be isolated. Since these interactions act on an (

) 2
Numerical simulations without exchange processes were conducted in Hilbert space by means of the Simpson program. [START_REF] Bak | A General Simulation Program for Solid-State NMR Spectroscopy[END_REF] The apparent exchange rate kT that includes contributions from both thermal activation and quantum tunneling can be expressed by: 23

𝑘 𝑇 = 𝐴 𝑅𝑇 ∫ exp ( -𝑊 𝑅𝑇 ) [1 + exp ( 𝑉 0 -𝑊 ℎ𝜈 𝑇 )] -1 𝑑𝑊 ∞ 0 , ( 3
)
where A is the classical Arrhenius-type pre-exponential factor, V0 the barrier height, and the tunneling frequency T is given by: 23

𝜈 𝑇 = 1 𝜋𝑎 √ 𝐸 0 2𝑚 , ( 4
)
where m is the total mass of the two tunneling protons and a the barrier width at height E0  V0. The integral in Eq.

Figure 1 .

 1 The parameters of the monoclinic unit cell are a = 8.92 Å, b = 7.83 Å and c = 9.43 Å, with  = 93.39° and a space group C2/c. 26 Four water molecules (shown in cyan) are found in each unit cell. Each water molecule is

Figure 1

 1 Figure 1 Crystal structure of Ba(ClO3)2•H2O, viewed down the c-axis. Water molecules are shown in cyan whereas all remaining barium, chlorine, and oxygen atoms are shown in orange, green and red, respectively. Water molecules are stacked along channels running down the c-axis, with alternating orientations in subsequent planes. Hydrogen bonds between water protons and oxygen atoms of chlorate anions are indicated by dashed lines. The lattice parameters obtained by neutron diffraction agree with the structure optimized with planewave-pseudopotential DFT methods, yielding a = 9.078, b = 8.061 and c = 9.591 Å, with  = 92.698°. The intramolecular 1 H-1 H dipolar coupling is ca. 30 kHz, whereas intermolecular couplings between neighboring H2O molecules are smaller than ca. 2 kHz.

Figure 2

 2 Figure 2 Numerical simulations of proton NMR spectra of an isolated H2O molecule in a polycrystalline Ba(ClO3)2•H2O powder spinning at a rate R = 10 kHz in a field B0 = 18.8 T (800 MHz), using intramolecular 1 H-1 H dipole-dipole and 1 H shielding tensors calculated by DFT. (a) Spectrum simulated assuming that the two proton shielding tensors are both orthogonal to the dipolar tensor and collinear with respect to each other, as expected when the CSA tensors are averaged by rapid exchange at high temperatures. (b) Spectrum that is expected at low temperatures where the dynamics are frozen, so that the CSA tensors of the two protons of each water molecule become "orientationally inequivalent", i.e., no longer collinear, although they retain the same principal components. One of the hitherto unknown "cross-term splittings" of ca. 3 kHz between the two sideband manifolds is emphasized by a tree diagram in red.

Figure 3 (

 3 Figure 3(Left) Experimental proton NMR spectra of polycrystalline Ba(ClO3)2•H2O acquired with a rotor-synchronized solid echo (R = 100 µs) while spinning about the magic angle at νR = 10 kHz and B0 = 18.8 T (800 MHz) over the temperature range 110  T  300 K. A recovery delay of 15 s between subsequent scans was used in all cases. (Right) Numerical simulations of spinning sideband patterns for water molecules with an intramolecular 1 H-1 H dipole-dipole interaction of 30 kHz and two orientationally inequivalent CSA tensors that are not collinear although they have the same principal components, with an isotropic chemical shift iso = 0 ppm, an anisotropy CS = -15.4 ppm, and an asymmetry ηCS = 0.12. The Euler angles relating the various interactions used in the simulations are given in Supporting Materials. The simulations were conducted in the 8 × 8 single-quantum Liouville subspace spanned by the basis {Ix, Iy, Sx, Sy, 2IxSz, 2IySz, 2IzSx, 2IzSy} relevant for the chemical shift and dipolar interactions acting on an initial state Ix + Sx. All simulated spectra result from powder averaging over 10 × 20 × 10 crystallite orientations. A Lorentzian line broadening of 900 Hz was applied to all calculated spectra.

Figure 4

 4 Figure 4Full widths at half height (FWHH) of the experimental MAS spectra of the left panel of Fig.3(black data points). Analogous FWHH of the simulations of the right panel of Fig.3(red data). Values were calculated with respect to the centerband. For T  140, the average FWHH for the two lines of the 'doublet' is reported.

Figure 5 (

 5 Figure 5(a) shows an Arrhenius plot ln(kex(T)) extracted from a comparison of experiments and simulations of Fig. 3 in the critical range 140  T  190 K represented by black

Figure 5 (

 5 Figure 5 (a) Temperature dependence of the exchange rates ln(kex) estimated by comparison of the experiments and simulations of Fig. 3 for 140  T  190 K (black dots with error bars). A linear fit is indicated by a black dashed line and yields Ea = 36.0  1.5 kJmol -1 and A = 8.4 × 10 15 (lnA = 36.7  1.0). Analogous data derived from the 2 H NMR study of Long et al. 19 are shown by cyan dots (Ea = 30.4 kJmol -1and A = 1.9 × 10 14 ). In grey, these latter data points are adapted for 1 H by assuming the usual kinetic isotope effect (Ea( 1 H) = Ea( 2 H)/2 1/2 = 21.7 kJmol -1 , keeping the same pre-exponential factor), in agreement with Ea( 1 H) = = 19.2 kJmol -1 estimated from the temperature dependence of T1( 1 H) by Silvidi.[START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF] The red line indicates an Arrhenius equation with a slope or thermal activation energy Ea = 44.2 kJmol -1 calculated by planewave pseudopotential DFT methods with an intercept or pre-exponential factor A = 1.4 × 10 18 chosen so as to match the kex of the 1 H case (grey data point) at high temperature (T = 300 K, left-most grey data point). The green line has the same slope as the red one, with a pre-exponential factor equal to that of Long et al.[START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] (b) Plot of the equation 3 (blue line) for the Bell-Limbach model where contributions of incoherent quantum tunneling are taken into account assuming our DFT parameters, and superimposed onto our experimental data points from (a) where six additional points have been added to indicate the temperatures for which simulated spectra are shown in Fig. 3. Data points are shown as empty dots without error bars in regions where the FWHH does not vary significantly with the temperature.

  Figure 5 (a) Temperature dependence of the exchange rates ln(kex) estimated by comparison of the experiments and simulations of Fig. 3 for 140  T  190 K (black dots with error bars). A linear fit is indicated by a black dashed line and yields Ea = 36.0  1.5 kJmol -1 and A = 8.4 × 10 15 (lnA = 36.7  1.0). Analogous data derived from the 2 H NMR study of Long et al. 19 are shown by cyan dots (Ea = 30.4 kJmol -1and A = 1.9 × 10 14 ). In grey, these latter data points are adapted for 1 H by assuming the usual kinetic isotope effect (Ea( 1 H) = Ea( 2 H)/2 1/2 = 21.7 kJmol -1 , keeping the same pre-exponential factor), in agreement with Ea( 1 H) = = 19.2 kJmol -1 estimated from the temperature dependence of T1( 1 H) by Silvidi.[START_REF] Silvidi | Protic spin-lattice relaxation in some hydrated barium salt crystals[END_REF] The red line indicates an Arrhenius equation with a slope or thermal activation energy Ea = 44.2 kJmol -1 calculated by planewave pseudopotential DFT methods with an intercept or pre-exponential factor A = 1.4 × 10 18 chosen so as to match the kex of the 1 H case (grey data point) at high temperature (T = 300 K, left-most grey data point). The green line has the same slope as the red one, with a pre-exponential factor equal to that of Long et al.[START_REF] Long | 2 H NMR Line Shapes and Spin-Lattice Relaxation in Ba(ClO3)2• 2 H2O[END_REF] (b) Plot of the equation 3 (blue line) for the Bell-Limbach model where contributions of incoherent quantum tunneling are taken into account assuming our DFT parameters, and superimposed onto our experimental data points from (a) where six additional points have been added to indicate the temperatures for which simulated spectra are shown in Fig. 3. Data points are shown as empty dots without error bars in regions where the FWHH does not vary significantly with the temperature.

  (a) shows the exchange rates predicted by the Bell-Limbach model in the case at hand when the height of the barrier is calculated by DFT. The same experimental data points of (a) are represented again in (b) by black dots, but six further six further points (open circles) have been added, corresponding to temperatures for which the spectra have been simulated in Fig. 3. These latter points do not show any error bars as they are related to motional regimes where the linewidth is not significantly influenced by changes in temperature, as shown in Fig. 4. The agreement between the Bell-Limbach model and our experimental data points is qualitatively satisfying over the entire temperature range. Figure 5(b) shows that a model that includes a thermally activated process at high temperatures and contributions of incoherent quantum tunneling at low temperatures describes our experimental data in the 140-190 K range in a satisfactory manner and agrees with our DFT calculations of the classical barrier height. For two interleaved spinning sideband manifolds to be observed in low-temperature spectra, as shown in Fig 2(b), two necessary conditions must be fulfilled: (i) the two protons must have two orientationally inequivalent (non-parallel) CSA tensors, (ii) they must have a mutual dipole-dipole coupling. The proof that these two conditions must both be fulfilled is given in the Supporting Information. Figure S3(a) shows a simulated MAS spectrum of a twospin system without exchange, like Fig 2(b). Figure S3(b) shows an analogous simulation where

Figure 6 Simulated 1 H

 61 Figure 6

  of water molecules confined in crystals of barium chlorate monohydrate over a wide temperature range 110  T  300 K. Numerical simulations that include proton exchange between two sites with non-equivalent (i.e., non-collinear) anisotropic chemical shift tensors reproduce the experimental spectral features remarkably well. The Arrhenius plot was compared with 2 H and 1 H NMR data of the same sample presented in the literature. The activation energy for the two-fold hop about the H-O-H bisector was evaluated by means of planewave-pseudopotential DFT methods. The Bell-Limbach model that includes contributions due to both a thermally activated processes and incoherent quantum tunneling is consistent with the experimental data over the entire explored temperature range. The non-Lorentzian lineshapes observed at low temperatures stem from cross terms due to interference between CSA and homonuclear dipolar interactions. Selective irradiation experiments prove that the two interleaved spinning sideband manifolds cannot be disentangled.

  initial state 𝝆(0) comprising only Ix and Sx, it is sufficient to consider an 8  8 single-quantum subspace spanned by the basis B = {Ix, Iy, Sx, Sy, 2IxSz, 2IySz, 2IzSx, 2IzSy}. In this basis, the corresponding superoperator Lex describing incoherent dynamics due to chemical exchange between the I and S spins is:

( 1 )

 1 was evaluated numerically over the interval 0  W  20 V0. The parameters utilized in Fig. 5(b) are A = 9.0  10 17 s -1 , E0 = 18.6 kJmol -1 , a = (/2)  1.61 Å and m = 2  1.673  10 -27 kg, corresponding to 2 proton masses traveling through a barrier with a thickness given by the semicircle spanned by a 180° rotation of an H2O molecule about the H-O-H bisector as taken from the DFT-optimized crystal structure. contract 'dilute para-water', grant agreement number 339754), and the French "Equipements d'Excellence" Paris-en-Resonance, contract ANR-10-EQPX-09.

  utilizing the Perdew-Burke-Ernzerhof (PBE) functional,[START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF] a cut-off energy of 80 Ry, 4  4  4 k-point grid, ultrasoft[START_REF] Vanderbilt | Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[END_REF] pseudopotentials for Ba and Cl atoms and norm-conserving[START_REF] Hamann | Conserving Pseudopotentials[END_REF] pseudopotentials for O and H atoms. After convergence, the geometry was characterized by total forces acting on all atoms smaller than 0.0038 Ry/Bohr. Cartesian coordinates and cell parameters are given in the Supporting Information. The optimized structure was subsequently used for a series of single-point energy (SPE) calculations to estimate the potential barrier of the 180° flip of one water molecule about the H-O-H bisector.

Acknowledgements

The authors thank Hans Limbach, Maurice Goldman, Bob Griffin, Jean Jeener and Malcolm Levitt for constructive comments and Nicolas Birlirakis and Mathieu Baudin for assistance with experiments. This work was supported by the CNRS, the European Research Council (ERC,


Ana Roxin

Elio Hbeich

Semantic interoperability between BIM and GISreview of existing standards and depiction of a novel approach

Keywords: Semantic interoperability between BIM BIM, GIS, Semantic Interoperability, Cyber-Physical Systems, Granularity

published or not. The documents may come

Semantic interoperability between BIM and GISreview of existing standards and depiction of a novel approach

Ana Roxin, Elio Hbeich

Introduction

Our today's society faces what we call the fourth industrial revolution (4IR) along with its impacts on our everyday lives. 4IR differs from the three previous revolutions because it not only addresses production automation but also knowledge automation. Sometimes referred at as "Industry 4.0" (the term was coined in Germany's manufacturing industries, a dozen of years ago), the epoch we are living is heavily impacted by the development of Cyber-Physical Systems (CPSs) in almost every activity sector. As the Web changed our lives 25 years ago, so the CPSs will also change the interactions we have with the physical world surrounding us. Usually defined as "physical and engineered systems whose operations are monitored, coordinated, controlled and integrated by a computing and communication core" (Rajkumar, 2010), CPSs can be applied in a multitude of application domains (e.g. agriculture, energy, buildings, manufacturing). Their design and implementation come with several issues that are usually best addressed through multi-disciplinary researches.

In this article, we consider the domain of urban processes' implementation, in the context of tomorrow's smart cities. Today's urban scopes usually come with a number of specific challenges that are too complex to address for a single stakeholder or territorial community. We are witnessing the occurrence of more and more individual behaviours and collective practices, along with innovative rules and norms (be it at a national or international level), all seeking to conceive multi-level and multi-scale solutions addressing the overwhelming challenges associated with the vision of smart cities. More specifically, when considering knowledge automation, smart cities become sandboxes for problem-solving, or incubators for intelligent approaches providing local answers to challenges such as climate change, energy efficiency or inhabitants' comfort. Cities and urban environments in general represent complex systems: they can hardly be reduced to their geographical scopes, and need to be interpreted from a multi-dimensional perspective encompassing their spatial, economic, social and cultural aspects. Complexity is pushed at an even higher level when taking into consideration the different standards and regulations that apply on each of the aspects listed above. Notably, regarding spatial aspects there are two main standard families that apply: a) Standards pertaining to Building Information Modelling (BIM), promoted by buildingSmart International (bSI) and the ISO TC 59 b) Standards pertaining to Geographic Information Systems (GIS), promoted by the Open Geospatial Consortium (OGC) and ISO TC 211. While GISs allow integrating different types of geographic information along with their attributes (e.g. raster images, digital maps), BIM aims at delivering methods for easing the management of building information thorough the built lifecycle (e.g. from design to demolition). While both standard families come with structured information models and processes for describing aspects from the considered domains, no links have been defined between the two worlds. Thus when it comes to implementing knowledge automation approaches in the context of smart cities, data must be seamlessly integrated into a system ensuring its consistent interpretation by the machine. In other words, interoperability must be reached among the models used for integrating the data. Or as such, a building's digital twin (as conceived with BIM standards) has no relation to a system's twin defining constraints and contexts the building has to integrate e.g. the urban scape. This prevents conceiving interoperable approaches for predictive maintenance, dynamic simulation or energy-efficiency improvement.

For addressing this issue, we present our approach for interoperability, based on meta-model federation. The article is divided as follows: section 2 introduces BIM and GIS information models as defined in the respective standard families, section 3 reviews existing standard approaches for interoperability, while section 4 describes BIM and GIS barriers. Our approach is discussed in section 5 and finally we conclude in section 6.

The need for interoperability

BIM information model

Building Information Modelling (BIM) is a 3D model-based process that gives architecture, engineering, and construction professionals the insight and tools to more efficiently plan, design, construct, and manage buildings and infrastructure. BIM model can be used for analysis to explore design options and to create visualizations that help stakeholders understand what the building will look like from start to finish. The model is then used to generate the design documentation for construction. Finally, BIM describes a method of work by which all relevant information for the life cycle of the building is integrated, administered and exchanged among the project participants. ISO 29481 (ISO 29481-1 2016) defines BIM as a shared digital representation of an object built to facilitate design, construction and operating process and form a reliable basis for decision-making. The first stage of BIM standardization was carried out in 1999 by IAI (now buildingSmart International) [START_REF] Eastman | BIM Handbook: A Guide to Building Information Modeling for Owners, Managers Designers[END_REF]. BIM relies on the following international standards:

 Information Delivery Manual (IDM) specifies how information is exchanged in a process. It is based on the ISO 29481 (ISO 29481-1:2016) standard and is defined as an interchange agreement. IDM is a natural language description of the exchange.  Model View Definition1 (MVD) describes the data model needed to meet the exchange requirements described in the IDM. The underlying methodology is described by Part 3 of ISO 29481 (ISO 29481-3:2010).  Industry Foundation Classes (IFC) [START_REF] Liebich | Industry Foundation Classes IFC4 official release[END_REF] represent the conceptual model for buildings and comprises all classes and relations for representing a building (ISO 16739 2013).

The IFC model is specified in EXPRESS and complies with ISO 10303 (ISO 10303-21:2002) also called STEP part 21 (STandard for the Exchange of Product model data). STEP focuses on the representation and exchange of product data and aims to integrate the processes of design, development, manufacture, and maintenance (see figure 1).

Figure 1: IFC layers of data schemas (ISO 16739-1) and modelled in the EXPRESS Schema (ISO 10303-11).

The GIS standard family

GIS allows capturing, storing, handling and analysing geographical data (Sahoo, 2017). The main international organization developing standards for geospatial information is ISO TC 211. ISO's Technical Committee 211 (ISO/TC 211) is dedicated to developing and deploying standards relating to geographic information. ISO/TC 211 specifies methods, tools, and services for data management, acquisition, processing, accessing, presenting, and transferring such data digitally (ISO 191xx series 2006).

The approach to conceptual modelling in the ISO 19100 series is based on the principles described in the ISO CSMF (Conceptual Schema Modeling Facilities) (ISO 191xx series 2006). This conceptual schema includes four levels: metamodels, conceptual (abstract) schemas, conceptual (applications) schemas and implementation schemas (see figure 2). The first contains the General Feature Model defined in ISO 19109, which specifies the concepts, terminology, operations, and assumptions needed to build the basic constructs in the Conceptual Schema layer level. The contents of the meta-meta model level is usually expressed in natural language and is not itself subject to standardization. Conceptual Schema layer contains the definitions of the concepts, terminology, operations and assumptions needed to construct application schemas. Application schemas define the types of features and processes that are instantiated to produce datasets of geographic information. Application schemas are expressed using syntax and semantics from one or more conceptual schemas. The "bottom" layer contains the actual data that is defined by the application schema at the application model level. 3. The need for interoperability

The concept of interoperability and its flavours

Defining "interoperability" isn't an easy task. Several definitions exist for this concept: the ISO alone holds more than a dozen standards, each coming with its own definition of "interoperability". The oldest definition of "interoperability" is from the ISO/IEC 1993 "Information Technology -Vocabulary -Part 1: Fundamental terms", and was updated in the vocabulary's 2015 version (ISO 2382(ISO :2015)). Interoperability is defined as the "capability to communicate, execute programs, or transfer data among various functional units in a manner that requires the user to have little or no knowledge of the unique characteristics of those units". This definition implies that when interoperable two systems can either exchange information or be accessed with a single method. The heterogeneity of definitions and interpretations harden the implementation of interoperable approaches in real-world applications and enterprises. In order to avoid potential ambiguities with "replace ability", "compatibility" is often used as a synonym of "interoperability" (ISO/IEC/IEEE 24765:2017, 3.2089). Given the above discussion about issues among GIS and BIM standard families, an approach for sustainable interoperability among those artefacts becomes more needed in related smart city knowledge automation applications. Indeed, as mentioned by the authors in (Dassiti et al. 2013), "in today's globally networked environment, one cannot achieve environmental, social/ethical or economic sustainability of any artefact […] without achieving ubiquitous ability of the artefact and its creators and users to exchange and understand shared information and if necessary perform processes on behalf of each other in other words, interoperate."

In order to further specify and tackle the interoperability issues among BIM and GIS, we follow the General System Theory (GST) abstraction (Von Bertalanffy et al.1969) and adapt it to the previous definition of CSPs. We thus consider BIM and GIS abstracted as systems comprising several parts, each part exhibiting some behaviour (that can be different from the overall system's behaviour). These behaviours and their related components, mechanisms and processes are monitored, managed and coordinated by some computer. Hence interoperability is achieved using standards that enable behaviours of parts of the system and the overall behaviour of the system to cooperate seamlessly in order to reach a common goal or function.

With these definitions and statements in mind, the next sections present existing levels of interoperability and discuss existing standard approaches for implementing interoperability.

Levels of interoperability

Existing standards identify three main levels of interoperability, namely: data, syntactic, and semantic interoperability. These layers are connected and build upon each other, lower levels providing elements required by upper levels functionalities [START_REF] Kubicek | Layers of Interoperability[END_REF]. The following figure illustrates those levels along with their definitions as pertaining to ISO standards.

Figure 3: Levels of interoperability

Sometimes referred at as physical interoperability, the issues pertaining to the data level of interoperability have been long resolved with the adoption of hardware standards such as Ethernet (IEEE 802.3,2018); along with standard protocols for lower layers of the ISO network architecture e.g. Coming back at our previous example, full semantic interoperability would allow having an OWL ontology preventing computer agents from deducing that "Dijon" is also the capital of "France". Following these definitions, semantic interoperability denotes the ability of applications and business partners to interpret exchanged data in a consistent way, implying explicit and formal structures. Such structures define the meaning of data elements and the relationship between them. As mentioned above, ontologies being "formal and explicit specifications of shared conceptualisations of a knowledge domain" (Studer et al, 1998) they represent the building blocks of semantic interoperability. Still, relying on ontologies doesn't lower semantic heterogeneity of the so conceived knowledge models. As an example, we can cite the numerous versions of OWL ontologies conceived for the IFC standard (Schevers & Drogemuller 2005), [START_REF] Beetz | IfcOWL: A case of transforming EXPRESS schemas into ontologies[END_REF], [START_REF] Zhang | Development of IFC-based construction industry ontology for information retrieval from IFC Models[END_REF], (Gao et al.2015). Indeed, while following the advice of Eastman et al [START_REF] Eastman | A Guide to Building Information Modeling for Owners[END_REF] suggesting that building data must be represented with Semantic Web technologies in order to reach semantic interoperability, all these ontologies have been defined independently from one another. No semantic links were defined to identify alignments between concepts and relations in those ontologies. Thus, no consistent interpretation can be delivered based on those ontologies solely. The need for defining links among existing knowledge models pertaining to BIM and GIS becomes urgent. And for doing so, the same approaches used for coupling models can be applied to ontologies. In this context, the ISO standard about the integration of industrial automation systems [START_REF]Industrial automation systems and integration -Concepts and rules for enterprise models. ISO 15704:2000 Industrial automation systems --Requirements for enterprise-reference architectures and methodologies SO 15745-3:2003 Industrial automation systems and integration --Open systems application integration framework --Part 3: Reference description for IEC 61158-based control systems CIBW78 -2019 -Submission 70[END_REF] defines three possibilities: models can be integrated, unified or federated. These approaches have been then withdrawn from ISO 14258 and integrated into ISO 15704 addressing the requirements for enterprisereference architectures and methodologies (ISO 15704, 2000). These three types of approaches were more recently considered as standard interoperability approaches in the context of ISO 11354, defining the Enterprise Interoperability Framework or EIF (ISO 11354, 2011). The sections below further discuss these three approaches, notably based on their specification in the EIF.

Standard Approaches for Semantic Interoperability

In the context of an integrated approach, all exchanged elements have to be represented with respect to a common form. Such common form must have an associated level of expressiveness allowing to capture the specific details of the elements exchanged, especially those impacting interoperability (ISO 11354, 2011). All elements and artefacts in the considered system or organization have to be described according to the common form, even if the latter isn't built upon an existing International standard. This approach is suitable when designing and implementing new systems rather than when reengineering existing systems for interoperability [START_REF] Métral | Ontology-based approached for improving the interoperability between 3D urban models[END_REF] ISO standards implementing model integration are the "Industrial automation systems and integration -Product data representation and exchange" (ISO 10303) and the ISO standard about "Enterprise integration -Constructs for enterprise modelling" (ISO 19440). Outside the ISO, model integration is also applied in the context of ebXML2 , a joint OASIS/UNCEFACT international initiative for enabling the consistent use of XML to "exchange electronic business data" and thus "facilitate open trade between organizations regardless of size". Unified approaches require a common meta-model. In its simplest version, such meta-model can be a reference vocabulary, while in a more advanced version it can represent a complete ontology. Defined as a meta-model, it allows establishing semantic equivalences among considered concepts or entities. All other considered models with their related syntaxes and semantics have to be mapped to the common meta-model. Using the common meta-model, a translation between the constituent models is possible even though they might encounter loss of some semantics or information. Unified approaches thus rely on model fusion (ISO 11354, 2011). Unified approaches are best suited for collaborative or networked environments, or in situations implying for a large enterprise to collaborate with several SMEs (ISO 11354, 2011). Examples of unified approaches are very common in the context of the ISO standard families from defined in the TC184 "Industrial data" (SC4) and "Interoperability, integration and architectures for enterprise systems and automation applications" (SC5) sub-committees. Most researches in the domain of interoperability also adopted unified approaches. For example, UEML (Unified Enterprise Modelling Language) aims at defining a neutral format at meta-model level to allow mapping between enterprise models and tools [START_REF] Métral | Ontology-based approached for improving the interoperability between 3D urban models[END_REF]. Federated approaches imply that no partner imposes their models, languages, or methods of work. Such approaches do not imply a common form or a common meta-model (ISO 11354, 2011). They mainly apply to contexts where the entities considered for interoperability rely on too different or too complex vocabularies or methodologies. In a federated approach, each entity needs to adapt its processes and methods. Computer agents are only provided with a priori information, about each entity along with their related capabilities. For reaching interoperability in such a context, mappings must be specified among input and output information of the considered entities or artefacts. Remaining inconsistencies must be manually addressed. Implementing successful federation among organizations or systems comes with more challenges than the two previous approaches. As an example of such implementation, we may cite the federation approach in ISO 16100 "Manufacturing Software Interoperability Services" (ISO 16100). The "Open systems application integration framework" and more specifically the profiles defined in its third part (ISO 15745-3:2003) e.g. process profiles, information exchange profiles and resource profiles, bring additional support for federation-based approaches.

Bringing semantic interoperability between BIM and GIS

While BIM comes with detailed 3D visualization and various functionalities to organize huge volumes of data related to buildings, GIS environments are highly customizable, well-equipped for multi-dimensional analysis, and ideal for projects involving multi-site environments. While BIM systems are best suited for managing data related to the buildings themselves, GIS applications pertain to the urban scope outside buildings. Even though one is usually struck by the differences among the methods and processes underlying both approaches, there is a general tendency of combining them in order to benefit from their cumulated advantages. Reaching a common vision in which BIM and GIS are complimentary to each other, would bring highly productive outcomes in the field of digital AECO (Architecture, Engineering, Construction and Operations). The integration of BIM and GIS can offer substantial benefits to manage planning processes during design and construction phases. While BIM systems focus on developing objects with maximum levels of detail for respective geometries, GIS are more focused on analysing the objects from the physical environment, based on different abstractions. Combining BIM and GIS processes and methods would allow a continuous and consistent interpretation of the data at different scales and from different point of views. However, reaching such vision comes with several challenges, the main ones being listed below: 1. Coordinate systems and spatial referencing: GIS use two dimensional real world coordinates (RWC 9), while BIM systems use three dimensional relative coordinates between objects, with a reference to RWC at root object. GIS is based on a global spatial reference systems and use boundary representation. BIM applications use local spatial reference systems. 2. Temporal aspects: In BIM applications, a building object is characterized by its geometrical representations and its geometrical and non-geometrical properties. Such object can have several geometrical representations, as they each correspond to a different point of view. Still, the BIM standards do not define any links between these geometrical representations and geometrical properties of the considered building object. Initially such permissiveness was wanted for BIM applications (in order to cope with how levels of detail are handled in GIS systems). But today, standards should restrict or specify explicitly the possible choices. The level of permissiveness allowed by today's standards hinders the efficient implementation of BIM ecosystems, as it all depends of the choices made at the level of software implementations 3. Semantics: BIM and GIS use different vocabularies to describe their entities and properties.

No equivalencies have been defined among these elements. While bSI and there is no define link between the IFC and GIS vocabulary has developed the bSDD (buildingSmart Data Dictionnary) listing all existing terms and properties in the IFC standard, there are no explicit links defined between the bSDD vocabulary and other similar initiatives such as the French standard XP P07-150 (AFNOR PPBIM), promoted in the context of CEN/TC 442 WG4. Such semantic links are essential for implementing consistent information exchanges based on the IFC format.

With respect to these challenges, next section will outline our approach for semantic interoperability between BIM and GIS ecosystems.

Our approach for achieving BIM/GIS semantic interoperability

Considering the above approaches, along with our application context e.g. knowledge automation in smart cities, approaches based on federation appear as the most suitable. Indeed, integrated approaches imply using one single common model according to which all other models are conceived and interpreted. As mentioned above (Table 1), these approaches are best suited when engineering novel CPSs, and fail in addressing all subtleties of existing CPSs. More specifically, in the context of our approach, two axis are considered for federation -horizontal, and vertical. For the first case, we consider relying on an existing approach namely the federated architecture for OWL ontologies or FOWLA [START_REF] Farias | FOWLA, A Federated Architecture for Ontologies[END_REF]. For the latter, Hobbs' granular partition theory [START_REF] Hobbs | Granularity[END_REF] gives several interesting perspectives and future work directions. Both approaches are discussed in the sections below. Following the database federation approach (Sheth et al, 1990), FOWLA is an approach relying on SWRL rules for federating autonomous ontologies (including TBox and ABox). The architecture contains two main components: The Federal Descriptor and the Federal Controller [START_REF] Farias | FOWLA, A Federated Architecture for Ontologies[END_REF].

The first is responsible of identifying missing concept instantiations and identifying new alignments (based on previously defined ones). The latter is mainly responsible of executing SPARQL queries. More specifically, it comes with a Rule Selector module that is responsible of selecting only the subset of SWRL rules that allow returning results pertaining to the considered SPARQL query. Granularity is the extent to which a system is composed of distinguishable pieces or grains. It can either refer to the extent to which a larger entity is subdivided, or the extent to which groups of smaller indistinguishable entities have joined together to become larger distinguishable entities. For example, a kilometre broken into centimetres has finer granularity than a kilometre broken into meters. Information granules, as the name itself stipulates, are collections of entities, usually originating at the numeric level, that are arranged together due to their similarity, functional adjacency, indistinguishability, coherency or alike [START_REF] Pedrycz | Granular clustering: A granular signature of data[END_REF]. The granular partition theory [START_REF] Hobbs | Granularity[END_REF] builds upon the classical extensional mereology, and considers that the world surrounding us can be represented through a global theory pertaining to First Order Logic theories. Granular approaches aim at extracting from this global theory, local theories that are less complex, easier to interpret and calculate. Thus, having P the ensemble of predicates available in a global theory T, D being the interpretation domain, for a local theory, one has to identify the ensemble of relevant predicates on P, namely R. This can be done by applying the indiscernibility relation defined by Hobbs [START_REF] Hobbs | Granularity[END_REF]: (∀x,y) [(x ∼ y)≡(∀p∈R)(p(x)≡p(y))] In order to best understand how this can be applied to our context, let us take an example. Consider planning a trip. In this case, the route one has to travel can be abstracted as a one dimensional curve. When considering an infrastructure use case involving for example works on the asphalt on the road, one can no longer approximate the road as a curve, but has to take into account its volume -it thus becomes a 3D volume. With the indiscernibility relation previously defined, one can identify predicates pertaining to the use cases considered. In the first one, two points in the asphalt, identified through their respective coordinates will be undiscernible. An example of a predicate pertaining in the context of this first use case would be the distance between one point on the road and the destination point. Granular computing is an approach orthogonal to existing modelling approaches. It allows separating one knowledge domain into smaller pieces of knowledge, by means of consistent and structured methods, thus building a granular perspective. This allows consistent reasoning on these smaller pieces of knowledge but also on the whole knowledge domain. Still, while several formal models of granularity have been defined in literature [START_REF] Mani | A theory of granularity and its application to problems of polysemy and underspecification of meaning[END_REF], [START_REF] Keet | A Formal Theory of Granularity -Towards enhancing biological and applied life sciencesinformation systems with granularity[END_REF], the different granular perspectives have to be explicitly and formally defined, with regard to the considered application domain. Moreover, in applications involving context awareness, one has to further study and specify the relation between knowledge granularity and context granularity. Given the above considerations, a first step in our approach addresses consistent semantic modelling of BIM and GIS information. Together with experts from the domain of BIM and GIS, the next steps of our work will investigate what alignments can be defined among BIM/GIS concepts and models (as defined in the respective ISO TC 211 and IFC ontologies). As such, the rules defined in the ISO 191xx standard family for application schemas (ISO 19109) and feature catalogues (ISO 19110) allow to represent IFC by means of UML. But as UML is not formal, additional alignments have to be investigated. More specifically, our future work will consider the following levels of alignments:

- Concepts from the IFC Shared Elements layer of data schemas have to be mapped to their respective equivalents in ISO 19130. With the above considerations in mind, future work to be done in the context of this approach also involves the following items:

-Missing ontologies: for example which ontology mediation will be used to establish compatibility on terminological level -Missing links: some can be identified fairly easily, others require exchanges with business experts and are more complex to define. -Granular approaches impact: the concept of granularity, seems intuitive and easy to implement, however the manner of ontologies, the associated levels and perspectives must be explicitly and formally specified by integrating the characteristics of the domain of knowledge concerned [START_REF] Livi | computational intelligence, and the analysis of non-geometric input spaces[END_REF]. In addition, when it comes to integrate granularity into application that handle business knowledge it is necessary to investigate, define and specify the granularity if knowledge and it is context.

Conclusion

In this article we aim at defining the interoperability issue among GIS and BIM systems, and specifying an approach addressing this issue. Our approach relies on Semantic Web technologies and granular approaches for performing two-axis federation. In our approach, we do not seek to merge BIM and GIS, neither to promote one over the other, hence we intend to reuse the FOWLA approach and its advantages in terms on lightly-coupled ontology federation. Granular approaches further help in conceiving and managing different abstractions of the same context or scape, which is highly pertaining to the urban environments considered by our application domain. The purpose of achieving interoperability between BIM and GIS is to specify and implement means to describe buildings along with their environment, at different scales.

Figure 2 :

 2 Figure 2: ISO TC 211 Conceptual model(ISO 191xx series 2006)

  Alignments among metamodels: the General Feature Model (GFM) of ISO 19109 has to be compared with the IFC elements contained in the core layer of data schemas of the IFC schema. IFC classes such as IfcKernel, IfcControlExtension, IfcProcessExtension, IfcProductExtension have to be mapped to their equivalents in ISO 19109 GFM. -Alignments among abstract conceptual GIS schemas and data schemas contained in the Ressource Definition layer of IFC: Several geometry and topology elements from the GIS temporal schema (ISO 19107) are equivalent to sub-classes of IfcDateTimeResource or IfcTopologyResource. Also several elements from ISO 19107 Temporal schema have equivalents in the IFC terminology notably subclasses of IfcGeometryResource and IfcPresentationAppearanceResource. IFC classes such as IfcGeometricConstraintResource or IfcGeometricModelResource have to be mapped to their equivalent concepts in ISO 19111 Geographic information -Spatial referencing by coordinates. -Alignments among application schemas in GIS and domain specific and shared elements IFC data schemas: IFC classes such as IfcKernel, IfcControlExtension, IfcProcessExtension, or IfcProductExtension have to be mapped to their equivalent concepts in ISO 19109 GFM.

  Semantic interoperability addresses the meaning of the messages exchanged and related issues have not yet been resolved by existing standards and approaches. Semantic Web standards and languages allow specifying such meaning, by means of formal and explicit specifications of conceptualisations e.g. ontologies. Considering the different Semantic Web languages existing, semantic interoperability comes with different flavours: Minimum semantic interoperability is enabled by the use of RDF ("RDF -Semantic Web Standards") and allows specifying the minimum knowledge than can be exchanged through a sentence e.g. what is expressed through the sentence itself. One simple example of such minimal semantic interoperability is a sentence (or message) specifying that the object "Paris" is linked to the object "France" through the property "is capital of". Such low level of semantic interoperability requires further manual and/or automated handling of the exchanged data.  Extended semantic interoperability allows defining a minimal ensemble of beliefs onto which two computer agents agree. Such ensemble of beliefs allows computer agents to make new deductions from the implicit facts contained in the message they exchange. Such level of semantic interoperability is enabled by the use of RDF Schema ("RDF Schema 1.1"). With RDFS, the elements forming a statement are identified by URIs (Unified Resource Identifiers), thus computer agents can dereference those URIs and access some shared RDFS-defined ontology specifying, for example, that a capital is a city, has a population, a name, etc. Extended semantic interoperability allows defining a common interpretation of the elements contained in the messages exchanged. There is no agreement upon what those elements may not be.  Full semantic interoperability is enabled by the usage of the OWL ontology language family.An OWL shared ontology can specify what computer agents may agree upon, while preventing them from making erroneous deductions. OWL allows specifying a knowledge conceptualization bounded to a given domain: the lower bound represents what computer agents are allowed to believe, while the upper bound identifies what they may not believe.

	Protocol).
	Data	Syntactic	Semantic
	interoperability	interoperability	interoperability
	"interoperability concerning the creation, meaning, computation, use, transfer, and exchange of data" (ISO/IEC 20944-1:2013, 3.21.12.4)	"ability of two or more systems or services to exchange structured information" (ISO 16678:2014, 2.1.18 )	"ability for data shared by systems to be understood at the level of fully defined domain concepts" (ISO/TS 18308:2004, 3.38)
		"interoperability such	"ability of two or more
		that the formats of the	systems or services to
		exchanged	automatically
		information can be	interpret and use
		understood by the	information that has
		participating systems"	been exchanged
		(ISO/IEC 19941:2017,	accurately" (ISO
		3.1.4)	16678:2014, 2.1.17)
			"interoperability so
			that the meaning of
			the data model within
			the context of a subject
			area is understood by
			the participating
			systems" (ISO/IEC
			19941:2017, 3.1.5 )

TCP/IP (RFC 791, RFC 793) and HTTP (RFC 2616). Syntactic interoperability addresses the syntax of messages exchanged among CSPs considered artefacts. The related issues have been resolves through the adoption of XML and related syntax standards e.g. HTML, WSDL (Web Service Language) and SOAP (Service Oriented Access

Table 1 :

 1 Comparison among the three standard approaches for semantic interoperability(ISO 11354, 2011)

		Integration	Unification	Federation
	Level of standardization	System level	Meta-level	Model level
		The form is not necessarily an	Establishes semantic	Links are explicitly and
	Advantages	international standard. All models are built and	equivalence allowing mappings between	formally defined at the level of the ontologies
		interpreted according to the	different models	themselves, and the sum

http://www.buildingsmart-tech.org/specifications/mvd-overview

http://www.ebxml.org/


Xiaoyu Wang
email: xiaoyu.wang@imt-atlantique.fr

Amine Bohi

Mariam Al Harrach

Mickael Dinomais

Julien Lefèvre

On early brain folding patterns using biomechanical growth modeling

Abnormal cortical folding patterns may be related to neurodevelopmental disorders such as lissencephaly and polymicrogyria. In this context, computational modeling is a powerful tool to provide a better understanding of the early brain folding process. Recent studies based on biomechanical modeling have shown that mechanical forces play a crucial role in the formation of cortical convolutions. However, the correlation between simulation results and biological facts, and the effect of physical parameters in these models remain unclear. In this paper, we propose a new brain longitudinal length growth model to improve brain model growth. In addition, we investigate the effect of the initial cortical thickness on folding patterns, quantifying the folds by the surface-based three-dimensional gyrification index and a spectral analysis of gyrification. The results tend to show that the use of such biomechanical models could highlight the links between neurodevelopmental diseases and physical parameters.

I. INTRODUCTION

Brain growth is accompanied by the folding of the cerebral cortex, which takes place during gestational weeks 16-40 [START_REF] Sun | Growth and folding of the mammalian cerebral cortex: from molecules to malformations[END_REF]. Recent studies have revealed that not only the molecular and cellular processes but also mechanical forces play an important role in the formation of the gyral and sulcal convolutions. For the latter, computational modeling is a key point to further understand brain folding mechanism. Based on a combination of two hypotheses for cortical folding, which are axonal tension and differential growth, very preterm neonatal brain growth is modeled by a morphogenetic growth of the cortical layer and a stretch-driven growth of the subcortical layer [START_REF] Budday | A mechanical model predicts morphological abnormalities in the developing human brain[END_REF]. The model and 2D simulation results show that mechanical stretch plays a crucial role in human brain development. Another human cerebral folding model (e.g., [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF], [START_REF] Tallinen | Gyrification from constrained cortical expansion[END_REF]), which is based on the hypothesis of differential tangential growth, can produce the characteristic patterns of convolutions on the brain surface over a developmental process. This model is guided by the use of 3D magnetic resonance images of a smooth fetal brain as an initial point. The 3D numerical simulations of the brain demonstrate that the relative tangential expansion of the cerebral cortex constrained by the white matter generates compressive stress, resulting in cusped sulci and smooth gyri similar to those in developing fetal brains.

However, the correspondence concerning the increase in brain volume, the development of the folding patterns and gestational age is not fully realistic in these models. To make the increase in the brain volume more in agreement with the biological measurement results (e.g., [START_REF] Chang | The assessment of normal fetal brain volume by 3-d ultrasound[END_REF], [START_REF] Hsu | Quantitative analysis of normal fetal brain volume and flow by three-dimensional power doppler ultrasound[END_REF], [START_REF] Armstrong | The ontogeny of human gyrification[END_REF], [START_REF] Gong | Fetal and fetal brain volume estimation in the third trimester of human pregnancy using gradient echo mr imaging[END_REF]), we propose a new brain longitudinal length (BLL) growth model, which is integrated into the human cortical folding model [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF], to improve the brain volume accuracy during the cerebral folding process.

The cortical folding patterns are influenced by various physical parameters, e.g., the initial cortical thickness [START_REF] Budday | A mechanical model predicts morphological abnormalities in the developing human brain[END_REF], the initial geometry [START_REF] Verner | A computational study of growth-driven folding patterns on shells, with application to the developing brain[END_REF] and the initial curvature of the surface [START_REF] Kroenke | How forces fold the cerebral cortex[END_REF]. A deeper understanding of brain physical parameters can significantly contribute to comprehend certain pathologies associated with abnormal cortical folding processes, which motivates us to study the effect of the initial cortical thickness on cortical folding patterns. In this study, we build numerical simulations of different initial cortical thicknesses using our combined model, which first allow us to visually remark the difference of the folding patterns. Then, we quantify these folding patterns by the surfacebased three-dimensional gyrification index and the curvaturebased spectral analysis [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (spangy) applied to adult brain size polymorphism[END_REF] to compare the complexity of the folding patterns.

II. METHODS

A. Biomechanical model of brain folding

In (e.g., [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF], [START_REF] Tallinen | Gyrification from constrained cortical expansion[END_REF]), the authors propose a human cortical folding model, which can mimic a realistic brain folding process. The brain is modeled as a neo-Hookean hyperelastic material solid. Brain growth is modeled by a relative tangential expansion of the cortical layer and the white matter layer, the cortical layer is assumed to grow more rapidly than the white matter layer. The brain model is discretized into highdensity tetrahedral finite elements for the simulations of the developing brain.

The initial condition for the growth of this brain model is derived from the 3D MRI of a 22 weeks' fetal brain, with the displacement and the velocity which are zero. The model obeys the Neumann boundary condition, in other words, the force is imposed on the brain model in finite element calculation. There are two main forces considered in the model. One is the elastic force, which is derived from the volumetric strain energy density of neo-Hookean and a deformation gradient. Another is the contact force, which takes place when a separation between a node and a triangle face at the brain surface is less than a threshold in order to prevent nodes from penetrating element faces.

The deformation gradient is defined in this model by F = A(G Â) -1 , which differs from the traditional definition of F = A Â-1 , where A is the deformed configuration of a tetrahedron, Â is the stress-free initial configuration of a tetrahedron. G is the relative tangential growth tensor, which associates with the distance of a tetrahedron from surface in material coordinates. Thus, when the relative tangential growth tensor is initialized, the brain model starts to grow and the deformation gradient is formed, the corresponding elastic force can be calculated. The contact force at the brain surface is obtained via penalty based vertex-triangle contact processing [START_REF] Ericson | Real-time collision detection[END_REF]. The resultant force is applied as the nodal force on each node of the mesh to produce the deformation of the brain model.

BLL growth model In order to explain the brain volume growth during the cortical folding process, Tallinen et al. propose a brain longitudinal length (BLL) growth model which is defined by BLL (mm) = 59/(1 -0.55t), where t parametrizes time and has a non-linear relation to gestational age (GA) as t = 6.926

× 10 -5 GA 3 - 0.00665GA 2 + 0.250GA -3.0189. t ∈ [0, 1] corresponds to GA ∈ [22weeks, adult].
According to the longitudinal length, the size of the initial brain will be scaled isotropically (in three dimensions). After measuring the brain volume of the simulation results, we find that this BLL growth model lacks brain volume growth during the folding process and cannot mimic a volume growth process consistent with the biological measurement results (e.g., [START_REF] Chang | The assessment of normal fetal brain volume by 3-d ultrasound[END_REF], [START_REF] Hsu | Quantitative analysis of normal fetal brain volume and flow by three-dimensional power doppler ultrasound[END_REF], [START_REF] Armstrong | The ontogeny of human gyrification[END_REF], [START_REF] Gong | Fetal and fetal brain volume estimation in the third trimester of human pregnancy using gradient echo mr imaging[END_REF]). Here, we propose a new BLL growth model. Referring to the normative fetal brain longitudinal length data in [START_REF] Kyriakopoulou | Normative biometry of the fetal brain using magnetic resonance imaging[END_REF], we find that a second-order polynomial is good enough to fit them and thus define the brain longitudinal length (BLL) as a function of gestational age (GA) from 22 to 40 weeks as

BLL (mm) = -0.067GA 2 + 7.16GA -66.05, (1)
with R 2 = 0.9987. The brain mesh is then scaled isotropically according to the longitudinal length of each gestational age.

In the cortical folding model (e.g., [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF], [START_REF] Tallinen | Gyrification from constrained cortical expansion[END_REF]), the relative tangential growth tensor is calculated for each tetrahedron by G = gI + (1 -g)n ⊗ n, which is perpendicular to the normal vector n of the tetrahedron. g is the expansion ratio of the gray matter relative to the white matter, which is defined by g = 1 + α 1+exp(10( y H -1)) , where α controls the magnitude of expansion and is defined as α = ( √ 8-1)t [START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF]. y is the distance from the top surface, which is calculated for four vertices of each tetrahedron and will be averaged. H is the cortical thickness. Considering that the cortical thickness changes little and increases steadily during the expanding process, we thus also propose to model the cortical thickness H as a linear function of time, as follows: H = H i +βt, where H i is the initial cortical thickness and β determines the magnitude of thickening. For a typical cortical thickness of 22 weeks' fetal brains, which is 2.5 mm, the corresponding deformed cortical thickness is approximately 2.8 mm in adulthood (t = 1), we thus define β = 0.3 in our simulations.

Then we combine the cortical folding model, the new BLL growth model and the model of the cortical thickness variation. Based on this combined model, we study in this work the effect of the initial cortical thickness on folding patterns.

B. Spectral analysis of gyrification

To quantify the cortical folding patterns of our simulations, we adopt the method for the spectral analysis of gyrification [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (spangy) applied to adult brain size polymorphism[END_REF]. This method performs a spectral decomposition of the mean curvature of the brain surface mesh based on the Laplace-Beltrami operator eigenfunctions. It allows to extend Fourier Analysis to more general domains such as graphs or surfaces and to produce power spectra. The Laplace-Beltrami operator is defined as ∆ M = div • ∇ M , where M is a Riemannian manifold. The eigenvalues of -∆ M are λ 0 = 0 ≤ λ 1 ≤ . . . and φ 0 , φ 1 , . . . are associated orthonormal basis of eigenfunctions, which satisfy

-∆ M φ i = λ i φ i .
(

) 2
The approach of the decomposition of the mean curvature (C) is described by

C = +∞ i=0 C i φ i , C i = M Cφ i , (3)
where C i is the Fourier coefficient of the curvature in the eigenfunctions basis φ i . A band power spectrum is defined by

BS C (0) = C 2 0 , BS C (k) = i k 2 i=i k 1 C 2 i , (4)
where k is the number of the band, B0 power corresponds to the first constant eigenfunction, other bands correspond to sums of powers in different frequency bands that defined via the wavelengths. i k 1 and i k 2 correspond to the numbers of eigenfunctions associated to the band interval limits.

The approach of the local dominant band segmentation, which is presented in [START_REF] Germanaud | Larger is twistier: spectral analysis of gyrification (spangy) applied to adult brain size polymorphism[END_REF], allows us to parcellate the cortical surface into 7 spatial frequency bands, from B0 to B6. The bands B1, B2 and B3 are related to global brain shape. The last three bands (B4, B5 and B6) are associated with foldrelated variations of curvature, which reveal the number of 1st, 2nd and 3rd order folds, respectively. The 3rd order folds locate where the absolute values of curvature are the biggest on the brain surfaces. In consideration that we want to quantify the folding variations induced by the difference in the initial cortical thickness, we focus thus on the folding power of B4, B5 and B6.

C. Surface-based three-dimensional gyrification index

In order to represent quantitatively the degree of cortical convolution, we use the surface-based three-dimensional gyrification index (3D GI), which is a global measure and defined as the ratio of the cortical surface area to the area of its smooth 'convex hull' (the minimum surface area needed to completely enclose the brain) [START_REF] Clouchoux | Quantitative in vivo mri measurement of cortical development in the fetus[END_REF]. For the simulations of different initial cortical thicknesses, we calculate and compare the 3D GI, which will be shown in the results section.

III. RESULTS

A. Validation of brain volume with biological time

Several numerical simulations of a 22 weeks fetal brain have been performed by using the combined biomechanical model described in Section II-A. We have then computed the fetal brain volume of the simulations with the initial cortical thicknesses varying from 1.48 to 2.98 mm. A comparison between the volume data of our BLL growth model, the volume data using Tallinen's BLL growth model and other validated data of the literatures (e.g., [START_REF] Chang | The assessment of normal fetal brain volume by 3-d ultrasound[END_REF], [START_REF] Hsu | Quantitative analysis of normal fetal brain volume and flow by three-dimensional power doppler ultrasound[END_REF], [START_REF] Armstrong | The ontogeny of human gyrification[END_REF], [START_REF] Gong | Fetal and fetal brain volume estimation in the third trimester of human pregnancy using gradient echo mr imaging[END_REF]) is shown in Fig. 1. Comparing with the volume data of Tallinen (cyan points), our volume data (blue points), fitted with a secondorder polynomial: BV (ml) = 0.59GA 2 -12.77GA+53.46, is much more consistent with other validated brain volume data of the literatures. At the same gestational age, the different points represent the brain volume of the simulations with different initial cortical thicknesses. The blue fitting curve of our volume data is consistent with these validated brain volume curves during gestational weeks, showing that our brain longitudinal length (BLL) growth model can mimic a realistic brain volume growth process.

B. Impact of cortical thickness on surface morphology

To understand the effect of the initial cortical thickness on the cortical folding patterns, we vary the initial cortical thickness in the combined model from 0.74 to 5.96 mm. The cortical thicknesses from 1.48 to 2.98 mm are defined according to normative fetal brain measurements. The other two cortical thicknesses (0.74 and 5.96 mm) are our hypotheses for fetal abnormal cortical thicknesses. The simulation results are shown in Fig. 2. At the 29th week of gestation, the brain surfaces already showed some clear linelike sulci. As the cortical thickness increases, the number of sulci decreases, the sulci became more isolated. Starting from 32 GW, smooth gyri begin to appear on cortical surface. With the increase in the cortical thickness, the gyri became wider and fewer. The qualitative analysis is consistent with previous works such as [START_REF] Budday | A mechanical model predicts morphological abnormalities in the developing human brain[END_REF] relating the cortical thickness and the folding patterns. A greatly increased cortical thickness associates with a severe decrease in number of gyri and sulci, which corresponds to the phenomenon of lissencephaly. The morphology with the initial cortical thickness of 5.96 mm in our simulations is similar to lissencephaly. On the contrary, polymicrogyria relates to an overly convoluted cortex, i.e., an increase in number and a decrease in size of gyri and sulci. The morphology with the thin initial cortical thickness of 0.74 mm resembles polymicrogyria. To analyze quantitatively the folding patterns, we first compute the 3D GI on the surfaces of the simulated brain (the initial cortical thicknesses are set between 1.48 and 2.98 mm), as is shown in Fig. 3. It can be seen that there is a clear increase of the 3D GI with the increase in gestational age, regardless of the cortical thickness. The difference in the 3D GI is small at the same gestational week. However, the effect of the cortical thickness on the 3D GI is irregular. The GI based on surface areas may not be able to describe the complexity of the sulci and gyri when the initial cortical thickness changes little. The 3D GI computed on cerebral hemispheres for 12 healthy fetuses [START_REF] Clouchoux | Quantitative in vivo mri measurement of cortical development in the fetus[END_REF] is plotted for comparison. We remark a good agreement on the 3D GI at early gestational ages, but bias appears at later gestational ages may be due to differences in individual brain development. In future work, we will use more fetal brain MRI data of different gestational ages and compare the simulation results with the data to further validate our model.

Secondly, for different initial cortical thicknesses (H i varying from 0.74 to 5.96 mm), we analyze the curvature spectra and compute the folding power of B4, B5 and B6 on the brain surfaces. The total folding power is shown in Table I. At the same gestational week, increasing the initial cortical thickness decreases the amount of growth in cortex which results in smaller total folding power, i.e., less total folds (1st, 2nd and 3rd order folds). Furthermore, we find that, for the thinner cortical thicknesses, the folding power of B6 accounts for an increasingly large proportion of the total folding power with the increasing gestational age. At each gestational week, the difference in B6 folding power is larger than that in B4 and B5 folding power between these cortical thicknesses, indicating that the small changes in cortical thickness have significant effect on the number of tertiary brain folds. These trends of B6 folding power for different cortical thicknesses are shown in Fig. 4. The thinner the initial cortical thickness, the larger the increase of the B6 folding power demonstrates.

IV. CONCLUSION

In this work, we have proposed to introduce a brain longitudinal length growth model derived from normative fetal brain measurement into a biomechanical cortical folding model. This has led to more realistic early brain development simulations validated by previous works of the literature. Furthermore, we have evaluated the impact of the initial cortical thickness onto the folding patterns through visual analysis and quantitative features. More specifically, experiments have shown that thinner cortical thickness leads to higher B6 folding power that corresponds to tertiary brain folds. Such biomechanical model may be used for future works to investigate causes of neurodevelopmental diseases associated with abnormal folding patterns of the cerebral cortex.

Fig. 1 .

 1 Fig. 1. The comparison of the brain volume development between our simulation results, simulations results of Tallinen's model and literature data.

Fig. 2 .

 2 Fig. 2. Cortical morphology with increasing gestational weeks for different initial cortical thicknesses.

Fig. 3 .

 3 Fig. 3. 3D gyrification index computed on each brain surface.

Fig. 4 .

 4 Fig. 4. Folding power of B6 for different cortical thicknesses.

TABLE I TOTAL

 I

	FOLDING POWER OF B4, B5 AND B6
		GW 22	GW 29	GW 32	GW 34
	H i = 0.74 mm	3.76	73.11	325.59	854.95
	H i = 1.48 mm	3.76	59.26	177.48	366.06
	H i = 1.98 mm	3.76	50.89	125.71	272.16
	H i = 2.48 mm	3.76	39.06	94.39	194.27
	H i = 2.98 mm	3.76	33.39	73.80	152.26
	H i = 5.96 mm	3.76	15.48	19.49	29.01

Aix Marseille Univ, CNRS, INT, Inst Neurosci Timone, Marseille, France.

[START_REF] Tallinen | On the growth and form of cortical convolutions[END_REF] LUNAM, Université d'Angers, LARIS -EA7315, F-49000, Angers, France.

ACKNOWLEDGMENT

The research leading to these results has been supported by the Fondation pour la Recherche Médicale (grant DIC20161236453).