Qin Qin

Xingyan Liu

Hanxiao Wang

Tingwei Sun

Fuqiang Chu

Lu Xie
email: xielu@ustb.edu.cn

Pascal Brault
email: pascal.brault@univ-orleans.fr

Qing Peng
email: andqing.peng@kfupm.edu.sa

Highly efficient desalination performance of carbon honeycomb based reverse osmosis membranes unveiled by molecular dynamics simulations

Keywords: Carbon honeycomb, Desalination, Molecular dynamics simulation, Temperature gradient

Seawater desalination is vital to our modern civilization. Here, we report that the carbon honeycomb (CHC) has an outstanding water permeability and salt rejection in the seawater desalination, as revealed by molecular dynamics simulations. More than 92% of ions are rejected by CHC at applied pressures ranging from 50 to 250 MPa. CHC has a perfect salt rejection at pressures below 150Mpa.On increasing the applied pressure up to 150 MPa, the salt rejection reduces only to 92%. Pressure, temperature and temperature gradient are noted to play a significant role in modulating the water flux.

The water flux increases with pressure and temperature. With the introduction of a temperature gradient of 3.5 K/nm, the seawater permeability increases by 33% as compared to room temperature. The water permeability of the carbon honeycomb is greater than other carbon materials and osmosis membranes including graphene (8.7 times) and graphyne (2.1 times). It indicates the significant potential of the carbon honeycomb for commercial application in water purification.

Introduction

The production of fresh water is a significant challenge for the increasing population in the society [START_REF] Werber | The Critical Need for Increased Selectivity, Not Increased Water Permeability, for Desalination Membranes[END_REF]. Seawater desalination is vital to sustain the modern civilization of human beings as one of the most important sources for obtaining fresh water from the oceans and brackish reservoirs [START_REF] Werber | Materials for next-generation desalination and water purification membranes[END_REF]. The reverse osmosis technology, which uses the pressure difference as a driving force to filter impurities, accounts for nearly half of the extent of water obtained from desalination [START_REF] Cohen-Tanugi | Nanoporous graphene as a water desalination membrane[END_REF], where the osmosis membranes play a crucial part, however, still suffer from instability and low efficiency.

The reverse osmosis membranes play an essential role in osmosis desalination and need to exhibit optimal stability, high water permeability and salt rejection for large-scale applications. Currently, commercial reverse osmosis membranes are widely used for salt rejection, however, the low water permeability (about 1 L/(m 2 •h•Bar)) and high capital costs have limited the development of such commercial membranes for water treatment [START_REF] Malaeb | Reverse osmosis technology for water treatment: State of the art review[END_REF]. Moreover, the conventional reverse osmosis membranes degrade in the presence of chlorine and are easily fouled by the organic matter [START_REF] Lee | A review of reverse osmosis membrane materials for desalination-Development to date and future potential[END_REF]. In this respect, many research efforts have been made to develop new nanoporous membrane materials with much higher stability and effective water desalination performance.

Carbon is known to be a flexible element with high stability. It can take form of various allotropes, such as carbon nanotubes (CNT), graphene and graphyne, with a varying arrangement of atoms. Recently, CNT [START_REF] Gethard | Water desalination using carbon-nanotubeenhanced membrane distillation[END_REF][START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF][START_REF] Thomas | Thermostat choice significantly influences water flow rates in molecular dynamics studies of carbon nanotubes[END_REF], graphyne [START_REF] Kou | Water permeation through singlelayer graphyne membrane[END_REF] and graphene [START_REF] Cohen-Tanugi | Water desalination across nanoporous graphene[END_REF][START_REF] Heiranian | Water desalination with a single-layer MoS2 nanopore[END_REF] have been extensively studied for water desalination using both experiments and simulations. In particular, the experimental studies have shown that the porous CNT membranes have a water permeability and salt rejection much higher than the currently used commercial membranes [START_REF] Gethard | Water desalination using carbon-nanotubeenhanced membrane distillation[END_REF], and the water flux increases with decreasing the membrane thickness [START_REF] Thomas | Thermostat choice significantly influences water flow rates in molecular dynamics studies of carbon nanotubes[END_REF]. Similarly, graphyne membranes have been widely investigated A c c e p t e d M a n u s c r i p t as porous carbon nanomaterial for desalination, and the results show that water can permeate swiftly through graphyne with high salt rejection [START_REF] Kou | Water permeation through singlelayer graphyne membrane[END_REF]. Recently, nanopores have been introduced into graphene's structure. Graphene can also be introduced into a variety of chemically functionalized nanopores (for example, MoS2 functionalized, hydroxylated and hydrogenated pores) and has been demonstrated to exhibit several orders of magnitude higher water permeability [START_REF] Cohen-Tanugi | Water desalination across nanoporous graphene[END_REF][START_REF] Heiranian | Water desalination with a single-layer MoS2 nanopore[END_REF][START_REF] Shi | Bilayer graphene with ripples for reverse osmosis desalination[END_REF][START_REF] Wang | Molecular dynamics study on water desalination through functionalized nanoporous graphene[END_REF]. With high water permeability and salt rejection, porous carbon nanomaterials have a high potential of commercial applications.

Pressure gradient is commonly used to generate controlled water flow in the nanoscale desalination membranes [START_REF] Werber | Materials for next-generation desalination and water purification membranes[END_REF]. In addition, the introduction of temperature gradient has been proved to induce water flow in the nanoscale porous nanomaterials.

Harvey et al. [START_REF] Zambrano | Thermophoretic Motion of Water Nanodroplets Confined inside Carbon Nanotubes[END_REF] studied the translation motion of the water nanodroplets driven by the temperature gradient in CNT. It was reported that the water flow imposed by the temperature gradient could be controlled [START_REF] Oyarzua | Carbon Nanotubes as Thermally Induced Water Pumps[END_REF]. The thermo-diffusion caused by the temperature gradient is known as the Soret effect or thermophoresis [START_REF] Dutrieux | On the Measurement of Positive Soret Coefficients[END_REF][START_REF] Duhr | Why Molecules Move along a Temperature Gradient[END_REF].

Carbon honeycomb (CHC), a new kind of nanoporous carbon nanomaterials, has been proposed by Park et al. [START_REF] Park | Electronic structure and mechanical stability of the graphitic honeycomb lattice[END_REF]. The special honeycomb-like structure of CHC has attracted a significant research attention since then. The sp 2 -sp 3 hybridized CHC was reported to have an outstanding stability with relatively small density and high specific surface per unit volume [START_REF] Hu | Threedimensional honeycomb carbon: Junction line distortion and novel emergent fermions[END_REF]. Krainyukova et al. [START_REF] Krainyukova | Carbon Honeycomb High Capacity Storage for Gaseous and Liquid Species[END_REF] experimentally developed stable CHC for the first time. Since then, the structure and properties of CHC have been widely explored [START_REF] Xie | Mechanical Properties of Vacancy Tuned Carbon Honeycomb[END_REF]. CHC has excellent performance in liquid storage [START_REF] Krainyukova | Carbon Honeycomb High Capacity Storage for Gaseous and Liquid Species[END_REF], along with enhancement of thermal conductivity [START_REF] Chen | Anisotropic thermal conductivity in carbon honeycomb[END_REF], electrical conductivity [START_REF] Wang | Semimetallic carbon honeycombs: new three-dimensional graphene allotropes with Dirac cones[END_REF] and mechanical properties [START_REF] Pang | Bottom-up Design of Three-Dimensional Carbon-Honeycomb with Superb Specific Strength and High Thermal Conductivity[END_REF][START_REF] Zhang | Nanomechanics of carbon honeycomb cellular structures[END_REF][START_REF] Qin | Anisotropic and temperature dependent mechanical properties of carbon honeycomb[END_REF]. For example, it was reported that the thermal conductivity of CHC in the cell axis direction was up to 100 W/mK with a density as low as 0.4 g/cm 3 , and the specific strength was higher than the other carbon materials [START_REF] Pang | Bottom-up Design of Three-Dimensional Carbon-Honeycomb with Superb Specific Strength and High Thermal Conductivity[END_REF].

A c c e p t e d M a n u s c r i p t

The aforementioned unique properties might be attributed to the distinctive structure: CHC essentially constitutes the rows of honeycomb holes. Rationally, it motivates the research community to further explore its properties and applications in desalination.

Herein, a CHC with sp 2 bonding in the wall and sp 3 bonding in the junction has been constructed, and the water desalination performance as a function of applied hydrostatic pressure and length have been investigated using the non-equilibrium molecular dynamics (NEMD) simulations. The effect of temperature and temperature gradient has also been investigated.

Models and Methods

Tilted view of CHC is shown in Figure 1(a). As a kind of asymmetrical armchair sp 2 -sp 3 bonded carbon honeycomb [START_REF] Wang | Atomistic Study of Mechanical Behaviors of Carbon Honeycombs[END_REF], it consists of graphene nanoribbons. The crosssection is a honeycomb-like structure based on a regular hexagon with side length of 5.8 Å. Each CHC unit consists of four pores (Figure 1(b)) with periodic boundary conditions in both x and y directions. The junction consists of an 8-carbon and two 5-carbon atom rings (as shown in Figure 1(c)). The water desalination system is composed of two graphene sheets as rigid piston at the two ends of the simulation system along the z-axis, a CHC, salt water in the feed side and fresh water in the permeate side (Figure 1 (d)). According to Robert S. Weatherup [START_REF] Weatherup | 2D Material Membranes for Operando Atmospheric Pressure Photoelectron Spectroscopy[END_REF], water will not escape from graphene and penetrate into the vacuum. The saltwater in the feed side consists of 12 Na ions, 12 Cl ions and 1200 water molecules (20×17×104 Å.). This, thus, produces a salt water area with salinity of 32 g • 𝐿 -1 , which is close to the salinity of seawater. Fresh water contains 300 water molecules with a size of 20×17×26 Å. Besides, period boundary conditions are used in all directions, along with two 35 Å vacuum layers placed at the two ends of the graphene piston so that no interaction occurs A c c e p t e d M a n u s c r i p t between the left and right pistons. Therefore, with the length of CHC as 5.8 nm, the size of the simulation box is about 20×17×258 Å. As the size of the simulation box becomes larger, it has a little effect on the obtained results. We studied the water flux by periodically replicating both x and y directions three times at 50 MPa. This way, the simulation box was enlarged by nine times. The difference in the water flux between the enlarged structure and pristine one is within 5%, implying that the result with such a simulation box is converged and independent of system size. In order to express the system structure more clearly, the simulation system with only 5.8 nm CHC has been shown in Figure 1d. The different lengths of CHC are shown in Supplementary Figure S1. NEMD simulations are employed to test the performance of CHC as water desalination membranes using the LAMMPS package [START_REF] Plimpton | Fast Parallel Algorithms for Short-Range Molecular Dynamics[END_REF][START_REF] Peng | Shockwave generates < 100 > dislocation loops in bcc iron[END_REF]. The post-processing including water flux, salt rejection and water molecules distribution is studied with

A c c e p t e d M a n u s c r i p t

OVITO [START_REF] Stukowski | Visualization and analysis of atomistic simulation data with OVITOthe Open Visualization Tool[END_REF][START_REF] Xie | Sensitive five-fold local symmetry to kinetic energy of depositing atoms in Cu-Zr thin film growth[END_REF]. We have used our resources to run three times for every case, and for each case, three MD simulations were performed by changing the initial conditions to achieve meaningful insights.

The interactions between the water molecules are described by the TIP4P model, and the salt ions are described by the non-polarizable model. The AIREBO potential [START_REF] Donald | A secondgeneration reactive empirical bond order (REBO) potential energy expression for hydrocarbons[END_REF] is utilized to describe the interactions between the carbon atoms. For other kinds of atoms, the Lennard Jones (LJ) and Coulombic terms are utilized [START_REF] Beu | Molecular dynamics simulations of ion transport through carbon nanotubes. I. Influence of geometry, ion specificity, and many-body interactions[END_REF][START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF][START_REF] Suk | Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations[END_REF]. The interaction energy () ij Ur is given as [START_REF] Frenkel | Understanding molecular simulation: from algorithms to applications[END_REF]:

( ) ( ) 12 6 ( ) 4 / / ij ij ij ij ij ij ij ij qq U r r r r      = - +   (1)
where ij  and ij  are the energy and length parameters, ij r is the distance between the atoms i and j, i q and j q are the charges on the two atoms and  is the dielectric constant.

The LJ parameters between the different elements are estimated by using the Lorentz Berthelot mixing rules. The LJ parameters are summarized in Supplementary Table S1.

The canonical (NVT) ensemble is used for the simulations with a Nose Hoover thermostat of 50 fs at T=300 K. The periodic boundary conditions are employed in all three directions of the simulated box, with a time step of 0.5 fs. Five different external pressures (ranging from 50 to 250 MPa) are applied on the graphene piston to achieve the water desalination process, and the pressure on the graphene piston on the fresh water side is keept at 0.1

MPa (1 atm). The simulation time is related to the applied pressure. The greater is the pressure, the shorter is the time taken by water to pass. 10,000,000 time steps are taken for each NEMD simulation.

It has been reported that the salt rejection is more important than water permeability for water desalination, which is calculated by [START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF]:

A c c e p t

Results and Discussion

Effect of pressure

Applied pressure has a significant influence on the performance of the water desalination performance of the membrane. Varying external pressures were applied to the two rigid graphene pistons in the cell axis direction so that the system is in a non- water flux is nearly linear with the applied pressure, which is consistent with the previous studies [START_REF] Heiranian | Water desalination with a single-layer MoS2 nanopore[END_REF].

Moreover, it has been reported that the salt rejection is more important than the water permeability for the desalination membranes [START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF]. The dashed line in Figure 2(b)

shows the salt rejection as a function of the applied pressure. The result shows that at low pressures, CHC has a perfect salt rejection with no salt ions passing before half of the water has flowed through CHC. On increasing the applied pressure up to 150 MPa, the salt rejection reduces only to 92%. This can be explained by considering the free energy profile of an ion and water molecule passing through the CHC membrane. As expected, the higher pressure provides enough energy to the salt ions to pass through the nanoporous A c c e p t e d M a n u s c r i p t membrane, however, at a very low extent. These results show that CHC could be an effective water desalination membrane alternative.

A few previous MD studies have shown that the water flux through CNT is affected by the length [START_REF] Thomas | Thermostat choice significantly influences water flow rates in molecular dynamics studies of carbon nanotubes[END_REF][START_REF] Su | Effect of nanochannel dimension on the transport of water molecules[END_REF]. In this study, the considered CHC lengths are 5.8, 11.6 and 58 nm respectively. The water flux of CHC with varying lengths as a function of applied pressure is shown in Figure 3(a). The decay of the water flux with increasing CHC length shows similar behaviour for different pressure values. As the length increases, the water flux is obviously noted to decrease. For CHC with a length of 58 nm, the water flux is noted to reduce by 17% as compared to CHC with a length of 5.8 nm at a pressure of 50

MPa. The reduction is observed to be higher on increasing the applied pressure. The maximum reduction of water flux is noted to be about 39% at 250 MPa. It has been reported that the use of the Nose-Hoover thermostats may overestimate the reduction in the water flux of the nanoporous carbon materials on increasing the length [START_REF] Thomas | Thermostat choice significantly influences water flow rates in molecular dynamics studies of carbon nanotubes[END_REF]. Therefore, the effective water permeability of CHC might be larger, if the similar effect occurs with CHCs.

A c c e p t e d M a n u s c r i p t At 300 K, the water permeability of CHCs with different lengths is described in Figure 3(c). As the CHC length increases, the water permeability decreases, exhibiting the same trend as reported for CNT for water transportation by Su [START_REF] Su | Effect of nanochannel dimension on the transport of water molecules[END_REF]. Remarkably, the water permeability of CHC with a length of 58 nm is 48% lower than that of 5.8 nm. It

shows the water flow through CHC is frictional [START_REF] Wan | Concerted orientation induced unidirectional water transport through nanochannels[END_REF].

The salt rejection as a function of the applied pressure is depicted in Figure 3(b).

Meanwhile, the membranes with CHC lengths of 11.6 nm and 58 nm length exhibit a perfect salt rejection, while the salt rejection of the CHC with 5.8 nm length is about 92% at 150 MPa. Similar salt rejection under other applied pressures is observed, indicating that the salt rejection is less affected by the CHC length.

Effect of temperature gradient

A c c e p t e d M a n u s c r i p t

Increasing temperature enhances the diffusion process of both water and salt ions.

It has been reported [START_REF] Oyarzua | Carbon Nanotubes as Thermally Induced Water Pumps[END_REF] that an axial temperature gradient can induce a fast and continuous water flow through CNT. Here, the effect of the temperature gradient on the performance of CHC for water desalination has been investigated. Simulations under the temperatures of 300 K and 320 K are taken as two reference groups. In order to discuss the effect of the temperature gradient, a different thermostat for each region is considered.

The atoms in the feed and permeate sides are heated to 320 K and 300 K, respectively.

The atoms in the middle region are controlled by the microcanonical (NVE) ensemble.

This process is noted to induce the temperature gradient of 3.5 K/nm in the water desalination system.

Water flux values at different temperature conditions are depicted in Figure 3 (a).

The results show that the larger water flux values are noted at 320 K than 300 K. The kinetic energy of the water molecules increases with temperature, which makes it easy for the water molecules to pass through. Thus, a higher water flux is noted at higher temperatures. The large pressure gradients can induce flow in CHC, so does the temperature gradient in this study. In order to optimally reflect the effect of temperature gradient on the water flux, we have further considered two temperature gradients (310-300 K and 330-300 K) at 50 MPa. The results show that the water flux increases on increasing the temperature gradient, as shown in Supplementary Figure S4.

In Figure 3(c), the water permeability of CHC has been plotted. The results show that the water permeability of CHC at 320 K is 13% higher than 300 K. The introduction of the temperature gradient improved the water permeability of CHC by 33%, indicating that the temperature gradient can promote the water permeability more effectively than simply raising the temperature.

A c c e p t e d M a n u s c r i p t

Previous study [START_REF] Insepov | Nanopumping Using Carbon Nanotubes[END_REF] has shown that the gases in CNT flow by applying traveling waves. Qiu et al. [START_REF] Qiu | Vibrating carbon nanotubes as water pumps[END_REF] found that the key to the water flow is the amplitude at the free end of CNT, which is induced by applying periodic force. The higher the amplitude, the greater is the water flux. Elton et al. [START_REF] Oyarzua | Carbon Nanotubes as Thermally Induced Water Pumps[END_REF] noted that the water flow in the temperature gradient direction within CNT can be caused by the difference in the amplitude of the oscillation between the high and low temperature zones. As the applied temperature gradient increases, the amplitude of the oscillations also increases, thus, causing the higher water flow rates. Here, the thermal Brownian ratchet concept [START_REF] Peskin | Cellular Motions and Thermal Fluctuations: The Brownian Ratchet[END_REF] is also suitable for the use of CHC for water desalination. By the induced temperature gradient, the water permeability also increases obviously. The introduction of the temperature gradient, thus, promotes the flow of water in CHC.

The Soret coefficient is calculated to quantify the effect of the temperature gradient on the water flux of CHC. The Soret coefficient is defined as [START_REF] Dutrieux | On the Measurement of Positive Soret Coefficients[END_REF]:

/ t S D D = (3)
where Dt and D are the thermal and isothermal diffusion coefficients respectively. The isothermal diffusion coefficient D is obtained from the classical mean square displacement:

2 1 1 lim ( ) (0) 6 N jj t j D r t r Nt → = =- (4)
where () rt is the displacement of the atom j at time t and N is the number of atoms. The coefficient of thermal diffusion is calculated from the equation:

t DT V =  (5)
A c c e p t e d M a n u s c r i p t

where V is the water flow rate determined by the water flux during the simulation. The isothermal diffusion coefficient is 2.8 nm 2 /ns, which is consistent with the previous reports for water under nanoscale confinement (0.9-8.8 nm 2 /ns) [START_REF] Zambrano | Thermophoretic Motion of Water Nanodroplets Confined inside Carbon Nanotubes[END_REF][START_REF] Striolo | The Mechanism of Water Diffusion in Narrow Carbon Nanotubes[END_REF]. The water flux induced by the temperature gradient under zero driving pressure is about 50 molecules/ns, corresponding to the water flow rate of 1.16 nm/s. Thus, the extracted thermal diffusion coefficient is 0.33 nm 2 /(ns K), which corresponds to a Soret coefficient of 0.12 K -1 . The observed value is one order of magnitude lower than the simulated coefficient for the confined water inside CNT [START_REF] Zambrano | Thermophoretic Motion of Water Nanodroplets Confined inside Carbon Nanotubes[END_REF], however, it is closer to the value for the binary Lennard-Jones liquids [START_REF] Reith | On the nature of thermal diffusion in binary Lennard-Jones liquids[END_REF] and the experimental results for the water-methanol mixtures [START_REF] Dutrieux | On the Measurement of Positive Soret Coefficients[END_REF].

The salt rejection as a function of both temperature and applied pressure is plotted in Figure 3(b). It is found to be 100 % for the applied pressure below 100 MPa. It is owing to the fact that the smaller values of pressure induce smaller forces on the salt ions to pass through CHC. The salt rejection makes a little difference at 150 MPa, where the temperature gradient shows a 100% salt rejection, which is higher than observed for 300 K and 320 K.

Porosity and water permeability

The water permeability is another important factor for the water desalination performance. Figure 4 shows the water permeability for various nanoporous carbon materials (hydroxylated and hydrogenated graphene [START_REF] Cohen-Tanugi | Water desalination across nanoporous graphene[END_REF], MoS2, graphene [START_REF] Heiranian | Water desalination with a single-layer MoS2 nanopore[END_REF], CNT [START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF] and graphyne [START_REF] Kou | Water permeation through singlelayer graphyne membrane[END_REF]) including CHC with L = 5.8 nm calculated in this study. CHC shows a better water permeability performance than the 2D membranes of graphene, MoS2, graphene with hydroxylated or hydrogenated pores and graphyne-5, with similar pore area. For 2D nanoporous materials, the pores tend to weaken the material by reducing its structural strength. It was reported that both strength and Young's modulus of the nanoporous graphene reduce on increasing the pore radius [START_REF] Cohen-Tanugi | Nanoporous graphene as a water desalination membrane[END_REF]. It is important for a A c c e p t e d M a n u s c r i p t material to be strong in order to maintain its mechanical structure under high hydraulic pressure inherent to the reverse osmosis desalination process. Thus, the porosity of such materials is generally small. This, thus, greatly limits water permeability of these materials. For instance, the water permeability of graphene with a 10% porosity is about 1000 L/(m 2 •h•Bar) [START_REF] Cohen-Tanugi | Water desalination across nanoporous graphene[END_REF].

Figure4. The performance of various membranes with similar single area for water desalination at 300 K. The single pore area of CHC is 78 Å2 with L = 5.8 nm. The pore area of the hydroxylated and hydrogenated graphene [START_REF] Cohen-Tanugi | Water desalination across nanoporous graphene[END_REF], MoS2 and graphene [START_REF] Heiranian | Water desalination with a single-layer MoS2 nanopore[END_REF] is 62, 50, 56 and 60 Å2 respectively, and the corresponding permeability is calculated assuming a 10% porosity. The pore area of CNT [START_REF] Shi | Bilayer graphene with ripples for reverse osmosis desalination[END_REF][START_REF] Shi | Bilayer graphene with ripples for reverse osmosis desalination[END_REF] [7] and Graphyne-5 [START_REF] Kou | Graphyne as the membrane for water desalination[END_REF] is 62 and 77 Å2. However, it is generally less critical to focus on the strength of CHC due to its unique nanoporous structure. The stress-strain curves along different directions show that CHC has outstanding structural strength (as shown in Supplementary Figure S5). In this study, CHC has a nanopore density of 1.15×10 14 cm -2 , while it is 1.7×10 13 cm -2 for graphene [START_REF] Cohen-Tanugi | Nanoporous graphene as a water desalination membrane[END_REF], 7.9×10 13 cm -2 for graphyne-5 [START_REF] Kou | Water permeation through singlelayer graphyne membrane[END_REF] and only 2.5×10 11 cm -2 for CNT [START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF] with similar pore area. For simplicity, the porosity is expressed by the ratio of the hole area and total cross-sectional area, where the hole area is calculated by the area of the inscribed circle of the hole. The calculated porosity of CHC is noted to be about 64%, whereas it is 43% for graphyne-5, 10% for graphene and only 0.5% for CNT.

Benefitting from the specific structure, the CHC membrane has a relatively high porosity, while maintaining a high stiffness. Therefore, CHC is noted to have a higher water permeability as compared to other materials. It has been reported that with a thickness of 3 nm and a diameter of 1.6 nm, CNT has a water permeability about 160 L/(m 2 •h•bar) [START_REF] Thomas | A computational assessment of the permeability and salt rejection of carbon nanotube membranes and their application to water desalination[END_REF], which is more than one order of magnitude lower than CHC with similar pore area.

Conclusions

A promising performance of CHC for water desalination has been predicted by means of molecular dynamics simulations. Benefitting from the specific structure of an array of honeycomb holes, the CHC membrane is noted to possess a relatively high porosity with high stiffness. To the authors' best knowledge, the water permeability of CHC is greater than other carbon materials and reverse osmosis membranes. The water transport through CHC is pressure dependent. Besides, the water permeability of CHC increased by inducing the temperature gradient, thus, providing a promising approach to improve the water flux. CHC shows high salt rejection for the applied pressure of less than 150 MPa. Overall, the newly developed CHC based water desalination membrane exhibits strong potential for seawater desalination and water transportation, along with application to the nanofluidic devices.

A c c e p t e d M a n u s c r i p t

Figure 1 .

 1 Figure 1. Simulation box and water desalination performance: (a) Angled and (b) vertical view of CHC; (c) Atom structure of the wall and junction of CHC; (d) Schematic of the simulation box consisting of a CHC, water, ions (Na + and Cl -) and two graphene sheets.

  equilibrium condition with a hydrostatic pressure difference across the CHC. The number of water molecules filtered through the CHC (length L = 5.8 nm) pores as a function of time for different pressure values is presented in Figure2(a). Each curve in Figure2(a)starts in the linear state with a near constant water flow rate. This suggests that the impact of the gradually increased salinity in the feed side could be negligible in this region. The water flow rate is estimated by the slope of each curve in Figure2(a). With time, the water molecules pass through CHC to the permeate side with increasing salinity in the feed side, following the osmotic gradient enhancement. Therefore, the driving force that causes the water molecules to flow is reduced, and the slope of the water molecules number vs. time curves decreases. At an externally applied pressure of 50 MPa, the water density and distribution of the oxygen atoms along the radial direction in a single pore are exhibited in Supplementary FigureS2 and S3, respectively. The results show that the water A c c e p t e d M a n u s c r i p t molecules are mainly distributed in the range of 1.6 to 2.4 Å from the center of the CHC pore.

Figure 2 .

 2 Figure 2. (a) The number of the water molecules filtered through the CHC pores as a function of simulation time at different pressures and (b)Water flux and salt rejection as a function of the applied pressure for CHC.

Figure 3 .

 3 Figure 3. Length and temperature dependence: Water flux (a) and salt rejection (b) as a function of the applied pressure for CHCs with different lengths at different temperature conditions; (c) Water permeability of CHC at different conditions.

A c c e p t e d M a n u s c r i p t

Acknowledgments

This work was supported by the National Key R&D Program of China (Grant No. 2020YFA0405700). Q. P. would like to acknowledge the support provided by the Deanship of Scientific Research (DSR) at King Fahd University of Petroleum & Minerals (KFUPM) through project No. DF201020.

The conflict of interest statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Tanguy Godquin
email: tanguy.godquin@ensicaen.fr

Morgan Barbier
email: morgan.barbier@ensicaen.fr

Chrystel Gaber
email: chrystel.gaber@orange.com

Jean-Luc Grimault
email: jeanluc.grimault@orange.com

Jean-Marie Le Bars
email: jean-marie.lebars@unicaen.fr

Applied graph theory to security: A qualitative placement of security solutions within IoT networks

Keywords: IoT, edge computing security, security as a service, centrality, dominating set

The rise of edge computing enables local network management. Services are no longer clustered inside the cloud but rather spread all over the whole network. In this paper, we propose a method for deploying security services within an IoT network according to devices capabilities. Our method models an IoT network as a weighted graph using device capabilities. Using the latter, we propose to identify the most suitable location for a security service using dominating sets and the graph weights. We present results obtained by the proposed method on an example of a smart city based on real data using network security functions such as an IPsec service. Results indicate an overall increase in the network security with minimal impact on the information flow while retaining reduced deployment costs.

Introduction

Since 2009, the estimated birth of the IoT by Cisco ISBG, the number of connected objects has been growing steadily, to reach in 2020, 25 (Gartner [1]), 50 (Cisco [START_REF] Evans | The internet of things: How the next evolution of the internet is changing everything[END_REF]) and even 200 billion devices according to Intel [START_REF]A Guide to the Internet of Things Infographic[END_REF]. The communication aspect of connected objects raises important security concerns. Often constrained in terms of hardware, software, and energy, a large majority of devices do not offer sufficient security to their communications. HP reports in 2014 that 70% of devices did not perform communication encryption [START_REF] Hp | HP Study Reveals 70 Percent of Internet of Things Devices Vulnerable to Attack[END_REF]. The Open Web Application Security Project (OWASP) foundation has compiled a top-10 list of mistakes to avoid during the conception of IoT applications, from password quality to the lack of hardware security [START_REF]OWASP Internet of Things Project -OWASP[END_REF]. Despite an increased use of security elements (mainly crypto-processors) and the ability of IoT devices to perform efficient cryptography [START_REF] Pereira | Performance Evaluation of Cryptographic Algorithms over IoT Platforms and Operating Systems[END_REF], a large number of devices are already deployed without these features. Faux et al. [START_REF] Faux | La sécurité à l'ère des objets connectés: Comment s'y prendre ?[END_REF] report that 84% of manufacturers do not evaluate the security of their products, which is even more alarming considering that 90% of IT security consultants are not trained to secure IoTs. These statistics confirm what Schneier states that "the market still largely rewards sacrificing security in favour of price and time-to-market" [START_REF] Schneier | Lessons From the Dyn DDoS Attack[END_REF]. The lack of security surrounding IoT devices is therefore the result of a design choice which leads to the deployment of highly vulnerable devices.

Regarding the Cloud, some concerns were expressed about the decentralization of cloud functionalities to the edge of the network. Edge computing paradigms relate to several computational paradigms such as fog computing or Multi-access Edge Computing (MEC). The fog computing is defined as follows: "a highly virtualized platform that provides compute, storage, and networking services between end devices and traditional Cloud Computing Data Centers, typically, but not exclusively located at the edge of network" [START_REF] Bonomi | Fog computing and its role in the internet of things[END_REF]. The Multi-access Edge Computing (MEC) refers to a paradigm initially focused on cellular networks that aims to "extend cloud computing capabilities to the edge of cellular networks" [START_REF] Porambage | Survey on Multi-Access Edge Computing for Internet of Things Realization[END_REF]. The telco aspect of MEC is particularly visible in the location of MEC nodes (or platforms). Unlike the fog computing, where the fog nodes are located at "any strategic location between end user device and cloud", the MEC platforms may be deployed at various locations within the Radio Access Networks (RANs) such as the Base Stations (which provide mobile services within a particular cell) [START_REF] Porambage | Survey on Multi-Access Edge Computing for Internet of Things Realization[END_REF]. The location of these platforms is critical when it comes to edge computing. As part of the MEC, telcos need to deploy services on specific RAN sites. Whereas in the context of fog computing, it will be essential to define what a strategic location is.

The migration from cloud to edge raises major concerns regarding the identification of suitable location for edge nodes and the allocation of services within them. The deployment of security services on every device is not realistic, either due to the deployment costs, the lack of capacities or restricted access to the device. Therefore, we ensure that every device from our network can directly communicate with a security node, a device that hosts a Network Security Function (NSF) such as for example IPsec services [START_REF] Battula | Network Security Function Virtualization(NSFV) towards Cloud computing with NFV Over Openflow infrastructure: Challenges and novel approaches[END_REF]. The aim of our work is to propose a method to deploy security services within an IoT network according to the capabilities of the device, including, but not limited to, processor, memory, storage, virtualization capacity and service-specific hardware requirements. The placement problem relates to the selection of nodes to host the security services. To limit the deployment cost of the services a minimal number of nodes have to be selected while providing the most suitable locations for the services. It forms an optimization problem with these two previous criteria. Furthermore, as stated above, our requirements dictate that each node on the network must be directly linked to at least one secure node. We formalize our problem through an undirected graph where V, the set of vertices, corresponds to the devices and E, the set of edges, refers to devices connections. Afterwards, each vertex is weighted to reflect the capabilities of the devices. The higher the weight, the more it has to be favoured, as it depicts a more compliant device, whereas nodes with a low weight will tend to be avoided. The requirement for direct access from each vertex of the graph to a security node implies that selected vertices forms a dominating set, i.e. a set where each vertex of the graph is either in the dominating set or directly linked to a node in this set. In addition, any dominating set, to solve our optimization problem, has to be of minimal size and favour heavier nodes. Our main contribution provides a greedy algorithm that identifies, from an undirected weighted graph, a dominating set satisfying our optimization criteria.

Several studies cover the service placement problem in IoT, edge or fog networks. Most of these studies focus on the placement of general services with numerous objectives. This may involve network concerns such as optimizing bandwidth and response times [START_REF] Jemaa | QoS-Aware VNF Placement Optimization in Edge-Central Carrier Cloud Architecture[END_REF][START_REF] Xia | Combining hardware nodes and software components ordering-based heuristics for optimizing the placement of distributed IoT applications in the fog[END_REF], the minimization of deployment costs [START_REF] Donassolo | Fog Based Framework for IoT Service Provisioning[END_REF], or the harmonization of performances [START_REF] Skarlat | Optimized IoT service placement in the fog[END_REF][START_REF] Skarlat | Towards QoS-Aware Fog Service Placement[END_REF]. Besides representing IoT networks in graph form, we have not found any approaches to the problem of service placement that rely on the use of graph theory algorithms and more specifically on graph centrality. While Banerjee et al. use graph centrality in IoT networks [START_REF] Banerjee | Optimality criterion for the insertion of multi-interface nodes to improve connectivity in heterogeneous IoT framework[END_REF], the authors do not perform a service placement and tackle a completely different problem known as vertex cover. Most related to our contribution is the work of Doriguzzi-Corin et al. [START_REF] Doriguzzi-Corin | Dynamic and Application-Aware Provisioning of Chained Virtual Security Network Functions[END_REF] which tackles the placement of chained security services using Integer Linear Programming (ILP) and proposes a heuristic for its resolution. This heuristic relies on a graph weighted according to network capacities as well as the computation of the shortest paths between its nodes with no use of graph centrality. The work from Doriguzzi-Corin et al. [START_REF] Doriguzzi-Corin | Dynamic and Application-Aware Provisioning of Chained Virtual Security Network Functions[END_REF] addresses the placement problem along a given path with fixed points of origin and destination. The objective is to secure these two points with the deployed service. In our case, we wish to secure the overall network, thus finding ourselves with multiple origin and destination points, which represents a different placement problem. Nevertheless, the method of the authors shows a very interesting approach, especially regarding the weighting used and the ILP formulation of their problem. We published in [START_REF] Godquin | Placement optimization of IoT security solutions for edge computing based on graph theory[END_REF] a partial solution of our placement IoT problem where the second criteria -we want the most suitable nodes -is partly considered. Indeed, we supposed that all devices had the same capabilities. Hence only the topology of the IoT network was used to define these nodes. In contrast, we consider in this paper a more realistic model integrating the capability of the devices. On the side of graph theory, we point how carefully the weighting must be performed, revealing a more challenging graph problem than we had before. A wrong choice can dramatically affect the solution. In this new study, we demonstrate that our methodology enables the deployment of security solutions on nodes suitable for hosting, thus drastically reducing deployment costs and more accurately reflecting the IoT environment while maintaining similar performance to previous work. While our security solution provides interesting benefits, other solutions can be considered, particularly regarding trust management such as [START_REF] Sato | Establishing Trust in the Emerging Era of IoT[END_REF] and [START_REF] Agostino Ardagna | Trustworthy IoT: An Evidence Collection Approach Based on Smart Contracts[END_REF].

This paper is organized as follows. Section 2 presents important contextual information regarding security constraints and attack threats on the network. A first approach is discussed, involving the application of edge computing concepts to IoT security. Section 3 provides an abstract methodology for qualitative positioning of edge computing security services on IoT networks. Section 4 outlines the proposed methodology behind a given security service and observes the benefits of this approach. Finally, we conclude this study in section 5 with a summary of our contributions supplemented by perspectives for further development of the model.

Security model

Despite an increased number of devices with on-board security features, the security is not always ensured. Features can be poorly implemented with outdated software or cryptographic algorithms. With a lifespan of roughly a decade once deployed, if not maintained correctly, a device will be subject to numerous vulnerabilities. A solution is to substitute all concerned objects by devices with suited capacities, however, it appears neither cost-effective nor realistic. A part of this problem have been addressed by Cloud Computing, in particular by centralizing services in powerful datacenters and servers. However, in the past few years, the trend has been to decentralize functions in the Cloud, a shift that is consistent with edge computing.

In this new paradigm, the gap between processing activities and information producers is reduced.

Among the most well-known examples of these mechanisms are the provision of video resources by streaming platforms on servers close to users who are most likely to consume such resources [START_REF] Bilal | Edge computing for interactive media and video streaming[END_REF] or processing tasks closer to information producers such as providing face recognition at the edge of the network [START_REF]Face recognition: Moving more and more to the edge -asmag[END_REF]. More generally, when it comes to fog computing, this frequently involves deploying Virtual Network Functions (VNF) to fog nodes. Regarding security, such functions (also called Virtual Network Security Functions: vNSF [START_REF] Battula | Network Security Function Virtualization(NSFV) towards Cloud computing with NFV Over Openflow infrastructure: Challenges and novel approaches[END_REF]) may be Virtual Deep Packet Inspector (vDPI), Virtual Firewall (vFirewall), Virtual Intrusion Detection System (vIDS) or even Virtual Encryption Proxy (vProxy) [START_REF] Farris | A Survey on Emerging SDN and NFV Security Mechanisms for IoT Systems[END_REF].

Threat model

It is acknowledged that no absolute security exists and the security of computer systems often consists in protecting a system against a specific threat. To propose a securing method, it is, therefore, necessary to define a model that reflects the capabilities and objectives of an attacker to protect a device effectively against that threat.

Stellios et al. [START_REF] Stellios | A Survey of IoT-enabled Cyberattacks: Assessing Attack Paths to Critical Infrastructures and Services[END_REF] discuss an attacker model suitable for the IoT composed of three categories: device accessibility, attacker capabilities and motivations. We consider an attacker based on the outsider model defined by Stellios et al., depicted on Figure 1.

The attacker has no physical nor proximal access to the object (in range of targeted device wireless protocols). As it is explained into the telecommunication standard ITU-T Y.4806, we can analyzed the security issues by their impact vector considering which layers is reachable by the attacker [START_REF]Security capabilities supporting safety of the Internet of things[END_REF]. In our working assumptions, the attacker can exclusively perform remote attacks. He can only access to his target from an intermediate device as shown by the arrows in Figure 1. An attack may pass through several nodes as it is the case for Attack Path 1 and Attack path 2. The attack can be prevented at any point along the path.

Resources of the attacker in the model are those of a regular person with moderate technical capacities (between neophyte and expert). Motivations of the attacker model depends on the targeted IoT network. We consider an occasional hacker or a cybercriminal whose resources remain moderate as per the definitions of motivations and attacker profiles provided by Mosenia and Jah [START_REF] Mosenia | A Comprehensive Study of Security of Internet-of[END_REF] and Onik et al. [START_REF] Onik | A Novel Approach for Network Attack Classification Based on Sequential Questions[END_REF]. Given a target device invisible from the Internet and the capacities of the attacker model, the considered threats use one or multiple intermediate devices to reach the target device, and the source is either on the Internet or in a local network. The presence of the attacker within the wireless range or his ability to physically access the target is not considered due to the higher cost of performing such attacks. Choosing this configuration is a trade-off on the costs of security in accordance with the proportion of attackers evicted.

Security as a Service Model for Edge Computing

To address this threat model, we propose to deploy NSFs at the edge of the network, thus benefiting from latency reduction, control over data processing and physical proximity of the objects. These security services can be deployed either as dedicated devices with sufficient capacities added to the network or as virtualized services in existing devices or gateways. Examples of such services are IPsec services, firewalls, intrusion detection systems, intrusion prevention systems [START_REF] Battula | Network Security Function Virtualization(NSFV) towards Cloud computing with NFV Over Openflow infrastructure: Challenges and novel approaches[END_REF][START_REF] Farris | A Survey on Emerging SDN and NFV Security Mechanisms for IoT Systems[END_REF]. Figure 2 illustrates the deployment of security solutions to prevent threats identified in Figure 1. It is noticeable that all three attack paths (red arrows in Figure 1) have been blocked as a result of the insertion of security solutions (green circles, D1, D2 and D4 in Figure 2) in specific devices on the network. There is no need to deploy security solutions on every device in the network, however, it must be located with regard to identified threats.

Model relevancy

To achieve the objective described in 2.2, the way security solutions are placed is of paramount concern. Indeed, it is necessary to maximize efficiency and coverage to minimize costs. Section 3 proposes a placement method to address the positioning of security solutions in an IoT network. Studying the most suitable locations for deploying security measures requires modelling an IoT network.

In the first part of our study, we proposed a simplified model for the placement of security functions. We intend to build upon this model to improve the quality of the deployment locations for security solutions. To improve our initial model, we identified two opportunities for improvements. One relating to the inherent capabilities of the devices and the other regarding the communication interfaces among them. The initial model assumed that no single device had sufficient security capabilities by itself and thus required to be improved to accommodate new security solutions. Similarly, we assumed that all the communication interfaces of our initial model had the same characteristics but a slight variation in them may have an impact on the transit of information.

To address the variation between device capabilities and communication interfaces, we propose to consider the inherent characteristics of the devices by introducing the concept of priority within communication interfaces.

Placement methodology

This section gives a general method to select which nodes will host the security services. Firstly, we introduce the notion of proximity with secure nodes and give the corresponding graph property. Then, we demonstrate that centrality metrics enable the vertices to be weighted considering the communications between the IoT nodes. Thus, we propose an algorithm that selects the vertices from a weighted graph by favouring the heavier nodes. Finally, we improve our method to consider a more realistic model where the characteristics of IoT nodes become part of the selection of security nodes.

Formalization of our problem

To provide a better security at the edge of the network, it is important to define the best positioning of this security.

Due to the myriad of connected devices, we propose to adopt a more abstract representation of IoT networks by transposing these networks into graphs with objects as vertices and communications as edges. This representation enables a higher degree of abstraction and takes advantage of graph theory concepts to IoT security.

In the following paragraph, d refers to the distance (measured in hop numbers) between an IoT device and its security service. Not being able to always provide security measures on objects (d = 0), which is the ideal solution, a compromise must be made. While the usage of the security capabilities of a direct neighbour (d = 1) may be acceptable assuming that an attacker does not interfere with this communication, proposing security methods at a higher distance (d ≥ 2) introduces security issues regarding the entities that are relaying the information. Moreover, a direct access (d = 1) follows the fog and MEC paradigms where the devices directly access the closest services [START_REF] Porambage | Survey on Multi-Access Edge Computing for Internet of Things Realization[END_REF].

Therefore, we seek to ensure that each object in our network is in direct contact with at least one device that offers security services that we call a security node. Let G(V, E) be an undirected graph where V and E correspond respectively to the vertices and edges of G, the neighbourhood of x ∈ V is defined as N [x]. The graph G will correspond to a IoT network where V will represent the devices and E will symbolize the connections between devices. Note that we will use the term node when we refer to the IoT while we use the term vertex of the graph to avoid the confusion between these data.

Modelling such conditions on the graph is equivalent to the identification of a dominating set. Indeed, a dominating set is a subset S of vertices of a graph G for which each vertex of G is either in S or adjacent to a vertex of S. To reduce the expense of converting devices into security nodes, one would like to identify a minimum dominating set. Finding a minimum dominating set is proved to be a NP-hard optimization problem where even the estimation of the minimum size remains difficult [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]. Aware of this situation, we do not attempt to identify the minimum dominating set; a minimal dominating set, i.e. a local minimum, is sufficient.

Minimal dominating set

We assume to have a vertex-weighting function. Algorithm 2 returns a minimal dominating set depending of these weights. The vertices are visited in ascending order and are kept out of the dominating set until there is no other choice but to place them into the dominating set, which has the effect of selecting heavier vertices. Doing so helps to select the heaviest nodes. Note that this is a general algorithm that can be used in other contexts; furthermore, its effectiveness strongly depends on the weight allocation.

We refer to F as the set of free vertices that have not yet been browsed by the algorithm. Once the vertex is scanned, it is either transferred to the dominating set (DS), to the set of covered vertices (C, vertices adjacent to vertices in DS) or to the set B corresponding to browsed vertices that will later be in C but require an adjacent vertex in DS. An adjacency list labelled N F (x) tracks the neighbourhood of the vertex x within F at any time. This suggested algorithm requires the following two functions (Algorithm 1) to be specified: F orcedDominant and P ropagation.

The F orcedDominant function specifies whether a vertex provided as a parameter should be added to a dominating set. To do so, it checks whether the vertex no longer has the possibility of being covered or whether one of its neighbours depends on it for its coverage. In this case, the addition into the DS set is required. The P ropagation function is performed when a vertex is added to the dominating set (DS). It covers all the neighbours of this vertex that have already been visited by the algorithm and that were on hold (B set). Move n from B to C

The proposed Algorithm (2) browses a given set of vertices, if the current vertex analysed is not covered yet but has at least one neighbour to be analysed (meaning it could be covered later by this neighbour), the latter is put on hold. When a vertex no longer has neighbours not analysed to provide coverage, the current vertex is added to the dominating set. Since the list of vertices is sorted by ascending weight given as input to the algorithm (vList), the dominating set discovered (DS) favours the heavier vertices. Lighter vertices are first placed in C or B which allows, later on, to integrate heavier vertices in DS.

IoT topology and centrality metrics

We proposed in [START_REF] Godquin | Placement optimization of IoT security solutions for edge computing based on graph theory[END_REF] a partial solution of our placement problem where all devices had the same capabilities and where only the graph topology was considered. This subsection outlines the method used. Note that in [START_REF] Godquin | Placement optimization of IoT security solutions for edge computing based on graph theory[END_REF] we were not talking about vertex-weight, only centrality value based on the graph topology.

Finding a small dominating set is not our only concern, we also wish to ensure that the positioning of the vertices making the dominating set is consistent with the transit of information in the graph which depends on its topology. Figure 3 illustrates a bad scenario where the vertices of a dominating set (in red) are not positioned along the information path.

Our first approach is to weight vertices according to their importance in the graph. The importance of a vertex can be obtained using centrality metrics. Depending on our security needs, some are more appropriate than others.

We limit our investigations on four of the most popular measures, namely degree centrality, eigenvector centrality, closeness centrality, and betweenness centrality which appear the more relevant for the IoT topology.

Degree centrality. The degree centrality is a measure of the importance of a vertex according to its degree.

Eigenvector centrality. The idea behind eigenvector centrality is that the importance of a vertex depends on the importance of its neighbours [START_REF] Ruhnau | Eigenvector-centrality-a node-centrality[END_REF][START_REF] Bonacich | Power and centrality: A family of measures[END_REF]. Given a graph G, the value of a vertex i ∈ G is calculated using the following Algorithm 2 Central Dominating Set Algorithm.

Input:

vList: ordered list of V Output:

DS: dominating set vertices

1: F ← V 2: DS, B, C, N F ← ∅ 3: while vList = ∅ do 4: a ← vList[0] 5: if a ∈ F then 6:
if ForcedDominant(a) then vList ← vList -a 18: Return DS equation where n is the number of vertices of G, A is the adjacency matrix of G and λ max is the largest eigenvalue (non-negative due to the Perron-Frobenius theorem):

C(i) = 1 λ max n j=1 A ij C(j).
Closeness centrality. The closeness value corresponds to the proximity of a vertex with all vertices of the graph [START_REF] Sabidussi | The centrality index of a graph[END_REF]. To compute this measure, all the shortest paths between the analysed vertex and the rest of the vertices of the graph are studied. The measurement of a vertex thus corresponds to the average distance between the vertex and the other vertices of the graph. The normalized mathematical representation of this measurement is the following where n is the number of vertices in the graph and g ij the geodesic distance between vertices i and j:

C (i) = n -1 n j=1 i =j g ij .
Closeness centrality is a graph metric adapted when one wishes to favour the vertices from which it is faster to reach all the other vertices of the graphs. It is used when positioning fire stations in cities [START_REF] Dehmer | Quantitative Graph Theory: Mathematical Foundations and Applications, Discrete Mathematics and Its Applications[END_REF] and could also be used for the positioning of update servers in IoT networks when required.

Betweenness centrality. This measurement corresponds to the number of times a vertex acts as a bridge on the shortest path between two other vertices in the graph [START_REF] Freeman | A Set of Measures of Centrality Based on[END_REF]. The value of a vertex i is expressed using the following equation where g jk is the number of shortest paths (or geodesics) connecting the vertices j and k, j = k, and g jk (i) is the number of these that pass through i:

C(i) = i =j i =k g jk (i) g jk .
Assuming the information in the graph flows along the shortest paths, the betweenness measures the contribution of a node to network traffic.

Depending on the centrality measure that is used, the importance of a vertex in a graph may differ. Figure 4 corroborates this statement by displaying the vertex with the highest value for the four defined measures on a given graph.

Real-world constraints consideration

Using centrality metrics exclusively is not enough to limit the deployment costs of security functions as much as possible while simultaneously protecting against the defined attacker model. To achieve this requirement, the devices suitable for the security functions must be identified. We define a priority value for nodes which includes several criteria such as processor capacity, memory, storage, bandwidth, specific hardware. We assume here that the vertices of our graph are weighted so that the most suitable nodes are the heaviest. We describe in next two subsections different approaches to compute the nodes weights according to their suitability.

The nodes with high priority value will then be prioritized for the selection of the dominating set. For our experiments and the evaluation of the dominating set, we introduce the concept of quality. The quality of a dominating set is the proportion of prioritized objects among this set. Using this methodology, nodes with sufficient hosting capabilities will be able to become security nodes through software updates, which is less expensive, while the others may be given a hardware update or even a complete replacement of the device.

To improve the quality of the selected dominating set, the computation of centrality values may include a priority value according to the characteristics of the devices and their communications. Node priority value will not be used directly to define the vertex weight function. We introduce an edge weighting function to induce the computation of centrality to consider priority node values and network topology.

Indeed, with our approach, edge weights impact the centrality value computation by favouring vertices which correspond to nodes with high priority value. As a result, the vertex weights are altered which implies the order of V in the Algorithm (2) to be changed. In next paragraphs we propose a method to compute the edge weight function from node priority values. Note that the priority value includes different criteria such as processor capacity, memory, storage, bandwidth, specific hardware.

While navigating the graph, if there are several paths available between two nodes, this weighting enables the selection of the path with the lowest weight. This approach consists in favouring communications between high-priority value nodes and thus increasing the quality of the dominating set identified.

Once a good dominating set has been identified, IoT devices corresponding to security nodes can either be updated to deliver desired security features or be replaced if their security capabilities are not sufficient.

Two devices interacting with each other are able to use those security nodes as gateways whenever it is necessary to benefit from the security functions offered by the latter ones. To ensure that communication has been carried out through one (or more if necessary) of these security gateways, the Proof Of Transit [START_REF] Brockners | Proof of Transit[END_REF] could be applied. For this procedure, one or more secrets are shared via Shamir Secret Sharing to nodes where one wishes to attest the passage of the communication. The security nodes would use their secrets (or sub-secrets) to sign communication packets flowing through them. Recipients could then verify that all previously imposed signatures had been applied to communications meaning that the security methods have been enforced.

Experiments over a specific use case

Dataset

Acquiring data related to communications between IoT objects is a difficult task. A large part of IoT datasets only contains values from various sensors.

Researchers working on the concept of Social IoT have made available a few datasets containing information on devices interactions [START_REF]Social Internet of Things[END_REF]. In our work, we use a graph generated from the OOR (Ownership Object Relationship) adjacency matrix containing information about public and private objects that we have filtered to keep only public static devices. The graph thus obtained is a non-oriented graph composed of 1,458 vertices and 35,657 edges. The vertices correspond to IoT objects while the edges are potential communications between these devices determined according to the available communication protocols of vertices (Bluetooth: 40 meters, WiFi: 400 meters and LoRa: 1,500 meters).

SIoT dataset is accompanied by a short description of the objects with the following information: the device identifier, the owner identifier, the category of the device, as well as information on the brand and model. These data are provided using simple identifiers ranging respectively from 1 to 16, 12 and 24. Details concerning the brand or model is not provided, however, categories are listed and included; Table 1 illustrates the distribution of these categories from the graph. In our methodology, Section 3, we discussed the need to perform a cleavage between devices according to their priority value. A distinction between two categories has been chosen to represent high-priority and low-priority value devices. In our study, depending on a required security service, an object will either have the ability to host that service or not. Introducing additional categories would allow variations in the hosting capability of a device. However, the difficulties covered in Section 4.7 and illustrated in Figure 7 would be increased exponentially depending on the number of categories.

The Smart Santander project, from which our dataset is derived, presents a smart-city use case for irrigation of parks and gardens where street lights are used as repeaters between the different communication interfaces of the network [START_REF] Smartsantander | [END_REF]; thus we have opted to define devices from the "Street Light" category as high-priority value nodes, hence prioritized. Moreover, street lights are no longer only considered as mere lights, many cities like San Diego, Chicago and Los Angeles invest in Smart-Street lights combining a variety of sensors [START_REF]Smart Streetlights Program -Sustainability -City of San Diego Official Website[END_REF][START_REF]Chicago smart lighting program[END_REF][START_REF] Garcetti | IES SALC Conference 2016 -Mayor Garcetti Speech[END_REF]. Now, they can host micro cell towers, provide WiFi coverage and even charge electric vehicles [START_REF] Garcetti | IES SALC Conference 2016 -Mayor Garcetti Speech[END_REF]. NEC offers smart street lights with embedded 5G base station and direct access to cloud applications [START_REF]Smart Street Lighting[END_REF]. Deploying security services on such devices would be relevant in a smart-city context. The ubiquity of street lights and their large number make them ideal targets to provide geographical coverage for different communication services, and therefore have the potential to host varied security functions.

Priority pre-processing

We fix two priority values L and H over the IoT nodes such that the smaller the value, the higher the priority of the node is. In our dataset the Street Light nodes will have high priority value H while the others have a low priority value L. For any vertex a ∈ V , we note P (a) its priority value. The edge weighting function is directly derived from P : for any e ∈ E, let a and b being its incidents vertices, W (e) = P (a) * P (b).

Table 2 displays the distribution of the weight of these categories over the considered graph. Weighting the edges of the graph to this extent affects the computation of the shortest paths used in centralities. The idea is to favour edges between nodes with high priority value while not overly penalizing the others. We show in our experiments that the weighting must be carefully chosen to offer the best results. The values of L were varied from 1 to 5 to identify which value gave the best results, these data are shown in Figure 7. We found that L = 1.2 and H = 1 are good choices. More details will be provided in Section 4.7.

Centralities

The effectiveness of an NSF relies on its location. While the NSFs perform operations on the network traffic such as filtering or encryption [START_REF] Battula | Network Security Function Virtualization(NSFV) towards Cloud computing with NFV Over Openflow infrastructure: Challenges and novel approaches[END_REF], it is clear that if such services are not located across the information flow, their efficiency will be impacted. Thus, the centrality metric that theoretically seems the most appropriate is the one that considers the flow of information in the graph, namely the betweenness centrality. Another advantage of betweenness centrality, also shared with closeness centrality, is to consider the weight of edges when computing shortest paths, meaning that a priority can be considered without any additional effort.

While computing shortest paths, edges with small weights (high priority value) will be preferred when navigating the graph. On the opposite, edges with a low priority value will be less likely to be selected as the paths will be longer. Using the weights in the centralities computation increases the importance (and therefore the value) of high-priority value nodes.

Dominating set selection

Once the values for all vertices in the graph have been calculated, our methodology suggests selecting a dominating set by favouring the addition of vertices with the highest priority values. More precisely, the selection operates as follows: the weakest value vertices are first put apart from the dominating set. A vertex is placed in the dominating set when it has no other options.

Evaluation

The algorithm introduced allows the analysis of the performance offered by applying various metrics. To reflect our requirements regarding the number of security nodes and vertices positioning, we use the following rating methods: the number of vertices in the dominating set, the quality of the dominating set, the protection provided through their positioning as well as a normalized value of protection independent of the graph size.

We defined the quality of a dominating set on Section 3.4 as the percentage of high-priority value nodes within the dominating set. The greater the percentage, the higher the overall quality will be considered. In our experiments, this refers to the proportion of devices from the "Street Light" category that are within the dominating set. Low-priority value nodes will either be upgraded or replaced with higher capacity versions to provide the requested service.

The protection provided by the vertex placement on dominating set is obtained by computing the percentage of vertex pairs from the graph that have at least one shortest path secured as follows: We define a shortest path to be secured if the vertices (sender and receiver) are adjacent or if each vertex of that pair is either member of the dominating set or adjacent to one. The illustration Figure 5 visually depicts this definition. Red vertices refer to dominating set vertices, while blue vertices are free to belong either to the set or not.

The normalized value of protection corresponds to the protection value normalized by the size of the graph over the size of its dominating set. It is computed using the following equation:

normalized value = protection * n |DS| .
This process is essential to consider the size of a dominating set when comparing the contributions of different methods of selection. Otherwise, a method returning a dominating set containing all vertices of a graph (100% protection ratio according to our evaluation) could be considered better than any dominating set.

Consider Figure 3, where we described the positioning of vertices not consistent with the information flow. According to our evaluation methodology, the protection provided by the dominating set placement on Figure 3 is 53%. Using the proposed approach for the placement of dominating set using the betweenness on the same graph, we obtain the graph in Figure 6 with a protection index of 100%. The dominating set discovered using our methodology is positioned across the information flow unlike in the one on Figure 3.

Results

According to the assumption that information flows along an optimal path, using betweenness appears more suitable for the placement of Network Security Functions. Algorithm 2 was used several times with lists of vertices ordered differently, the resulting dominating sets were then examined. We focused our experiments on the four most popular metrics, namely closeness (CC), eigenvector (EC), degree (DC) and betweenness (BC). The proposed algorithm favours selection of vertices located at the end of the list provided, therefore these listings were sorted by increasing order of centrality. The dominating sets obtained are analysed according to our evaluation criteria (size, quality, protection provided by the set and protection brought by one vertex) to be compared with those identified from randomly ordered lists (R) of vertices. Comparison to randomly ordered lists of vertices highlights the importance of vertex order in the selection of dominating sets. The results of this comparison are summarized in Table 3 below while Figure 8 summarizes for each category of devices, their proportion among the dominating sets. The effect of weight ratio has been analysed using betweenness (BC) and closeness (CC) centralities. The results are plotted in Figure 7 where each curve corresponds to an evaluation criterion according to the ratio between weights of low (L) and high (H) priority nodes: L/H.

Analysis on IoT dataset

Results from this study indicate that the betweenness and closeness centralities provide smaller dominating sets. It is noticeable that using centrality values in general offers better performance overall. Table 3 presents only two random tests by choice of readability of the data. The dominating sets identified using a collection of random lists of vertices benefit from similar performance to those depicted in the last two columns of Table 3 ranging from 180 to 210 vertices with a proportion of protection from 81% to 99.91%. Dominating set sizes are about twice as small using betweenness and closeness values as when vertices are sorted in random order. The gain observed is not limited to the number of vertices of the dominating sets but also their positioning. While the majority of the results in Table 3 reveals protection level approaching 100% (only reached with betweenness on weightless graph), it should be remembered that this value depends on the number of vertices in the dominating set. A larger dominating set is more likely to have a high protection value. Normalizing the measure according to the number of vertices of sets is therefore important for an accurate representation of gains. Protection index provided by a vertex from our IoT reference weightless graph using betweenness is almost twice as high as for random lists. The metrics of betweenness and closeness centralities stand out particularly from this analysis by offering a gain of one-third superior to the degree (highest index of all other results).

The dominating sets identified by betweenness and closeness over the weightless dataset, although highly similar from a performance point of view, are not identical. There are 18 vertices that are not shared between them, this similitude lies in the calculation methods of these metrics using both shortest paths within a graph.

Observations on the non-weighted graph can be reported over the weighted graph. As expected, there is a slight increase in size of the dominating set, due to the additional constraints added on the selection process. With the addition of the weighting, the protection value is reduced. It means that fewer shortest paths in the graph are considered secured as depicted in Figure 5. In our dataset consisting of 1,458 nodes, 35,657 edges and 1,062,153 shortest paths (potential communications), only 2,275 communications will have additional cost in their communications.

Increasing the quality of the dominating sets identified with weighting is particularly interesting when reducing the cost to deploy security functions. We defined the quality as the proportion of devices that can directly host the desired security function. Software deployment may be sufficient on this category of devices. Assuming that software deployment is more cost-efficient than upgrading or replacing devices, increasing the quality of the dominating set can significantly reduce the deployment costs of a security function.

It would not be possible to achieve a quality of 100% in our dataset since all 506 street lights nodes do not form a dominating set. Yet, weighting the graph improves the quality of the dominating sets while limiting the number of nodes added and additional communication costs.

Figure 7 shows the impact of the ratio between the low priority value and the high priority value over the number of vertices in the dominating set and the quality of the solution.

Closeness centrality attaches increasing importance to vertices that can spread information the fastest in a graph while betweenness reveals vertices that are located across the information flow. We are more likely to select betweenness to place NSFs across the information flow. Closeness cannot be excluded, however, it could very well correspond to a use case requiring the propagation of information through an IoT network.

We have decided to distinguish the measurement choice from the proposed algorithm and our methodology in order to allow the user to decide which option is best suited to his security needs.

Conclusion and future work

We propose a strategy to effectively deploy security solutions within an IoT network while minimizing costs. With the help of graph modelling, we reformulate the security constraints as the identification of dominating set where we favour the heavier vertices. Using centralities to compute the dominating set greatly reduces its size while accepting some adjustments using the graph weighting. We modelled the compliancy (ability to host a requested service) by putting a priority value on the nodes. Our solution has a minimum impact on the information transit while respecting the assumption that it flows optimally through the graph. We demonstrate that classifying nodes into two categories, high and low value, requires a thorough choice of the values used. Increasing the number of categories could offer better results. However the weighting will be exponentially more challenging. Introducing weights in the selection of the dominating set significantly increases its quality while preserving the low impact on communications and limiting the number of nodes added. Adding weighted edges on the graph therefore reduces the cost to deploy security functions by favouring devices that already have sufficient hosting capabilities.

Interesting perspectives could be to improve the proposed model by including other IoT features (more detailed device capabilities, risks or even bandwidth) to achieve a more realistic model, more specifically regarding weight assignment on the graph.

Figure 1 :

 1 Figure 1: Different examples of attack paths in the IoT context.

Figure 2 :

 2 Figure 2: Example of applied security solutions (green plain nodes) over our considered model.

Figure 3 :

 3 Figure 3: Dominating set (red vertices) located off the data flow.

Algorithm 1 Functions used on algorithm 2 1: 3 :T rue 4 : 6 :T rue 7 : 9 :

 1234679 function ForcedDominant(node) 2: if N F (node) = ∅ then Return for all n ∈ N [node] ∩ B do 5: if |N F (n)| = 1 then Return Return F alse 8: function Propagation(node) for all n ∈ N [node] ∩ B do 10:

7 :Propagation(a) 9 : 11 : 12 :

 791112 Move a from F to DS 8:else if |N F (a)| = 1 then 10: b ← N F (a)[0] Move a from F to C Move b from F to DS

Figure 4 :

 4 Figure 4: Difference of centrality values on a given graph.

  , b) = 1 if there is a secured shortest path between vertices a and b S(a, b) = 0 otherwise.

Figure 5 :

 5 Figure 5: Illustration of secured paths

Figure 6 :

 6 Figure 6: Dominating set (red vertices) identified using our methodology on Figure 3 graph.

Figure 7 :

 7 Figure 7: The impact of weight ratio on the evaluation criteria of the dominating set between prioritized and non-prioritized nodes.

  (a) Betweenness (114 nodes) (b) Closeness (132 nodes) (c) Degree (155 nodes) (d) Eigenvector (176 nodes) (e) Random (seed=82) (191 nodes) (f) Random (seed=45) (212 nodes)

Figure 8 :

 8 Figure 8: Quality of central dominating sets over our weighted dataset (weight ratio: 1.2).

Table 1 :

 1 Distribution of device categories in our graph

	Category	Number of devices Proportion
	Point of interest	95	6.5%
	Environment and weather	140	9.6%
	Indicator	10	0.7%
	Street Light	506	34.7%
	Parking	677	46.43%
	Alarms	30	2%

Table 2 :

 2 Weight distribution of edges in the graph

	Weight Number of edges Proportion
	1	961	2.7%
	1.2	9,088	25.49%
	1.44	25,608	71.82%

Table 3 :

 3 Comparison of central dominating set results between ordered and randomly ordered vertex lists.

				Centralities		Random order
	Graph	Metric	BC	CC	DC	EC	R seed=82	R seed=45
	IoT dataset	DS size	107 (7.33%)	107 (7.33%)	155 (10.63%)	176 (12.07%)	191 (13.10%)	212 (14.54%)
	weightless (1,458	Protection (%)	100	99.92	99.95	99.91	90.87	99.16
	vertices)
		Normalized	13.63	13.62	9.402	8.277	6.937	6.820
		protection
		Quality	47.7% 48.6% 44.5% 46.0% 51.3%	52.4%
	IoT dataset	DS size	114 (7.819%)	132 (9.053%)	155 (10.63%)	176 (12.07%)	191 (13.10%)	212 (14.54%)
	weighted (1,458 vertices)	Protection (%)	99.79	99.79	99.95	99.91	90.87	99.16
	weight
	ratio:	Normalized	12.76	11.02	9.402	8.277	6.937	6.820
	1.2	protection
		Quality	71.9% 81.8% 44.5% 46.0% 51.3%	52.4%

Acknowledgment

The authors would like to thank Paul Dorbec and Nancy Perrot for their help regarding dominating sets and their constructive remarks.