Arsenia Chorti
email: arsenia.chorti@ensea.fr

Comparison of Short Blocklength Slepian-Wolf Coding for Key Reconciliation

In this work, we focus on Slepian-Wolf (SW) coding in the short blocklength for reconciliation in secret key generation and physical unclonable functions (PUF). In the problem formulation, two legitimate parties wish to generate a common secret key from a noisy observation of a common random source in the presence of a passive eavesdropper. We consider three different families of codes for key reconciliation. All the selected codes show promising performances in information transmission in the short blocklength regime. We implement and compare the performance of the codes for SW reconciliation in the terms of reliability and decoding complexity for different blocklengths.

I. INTRODUCTION

Wireless communication technologies have become an essential part of our everyday lives. As more and more data are being transmitted over wireless channels, besides the reliability of the transmitted information, guaranteeing the security of information transferred has become a challenging issue. Physical layer security (PLS) has been investigated in the recent years to provide the security in information exchanged information by exploiting the properties of the physical layer [START_REF] Chorti | Physical layer security: A paradigm shift in data confidentiality[END_REF]- [START_REF] Chorti | Context-aware security for 6G wireless the role of physical layer security[END_REF].

In the context of PLS, in this paper we consider the problem of coding for secret key generation (SKG) from a common source of randomness observed by two communication parties 1 . In the problem formulation, the observation of two legitimate parties, who wish agree on a common secret key, differ due to noise or other factors. The central idea of reconciliation coding for SKG is that the parties should be able to correct the mismatches, by sending some helper data trough the public communication channel. Of course, the helper data can also help an eavesdropper (Eve) to have a better guess about the respective measurements.

Therefore, in coding for SKG, the target is to design a coding scheme which requires the minimum amount of helper data (information leakage), while guaranteeing that the legitimate parties will be able to reconcile the mismatches with a high probability (reliability). The coding for key reconciliation is also referred to as Slepian-Wolf (SW) coding in the literature. Motivated by the increasing importance of low latency communications in recent years, we focus on the short blocklength SW coding in this contribution. 1 The proposed codes can be further applied in the reconciliation step in the context of physical unclonable functions [START_REF] Mitev | Authenticated secret key generation in delay-constrained wireless systems[END_REF].

The problem of coding for SW reconciliation has been studied in the literature for long blocklengths. The SKG based on polar codes have been studied in [START_REF] Chen | A robust SRAM-PUF key generation scheme based on polar codes[END_REF], [START_REF] Hentilä | On polar coding for finite blocklength secret key generation over wireless channels[END_REF]. In [START_REF] Maes | PUFKY: A fully functional puf-based cryptographic key generator[END_REF], [START_REF] Forchhammer | Rate-adaptive BCH coding for Slepian-Wolf coding of highly correlated sources[END_REF], Bose Chaudhuri Hocquenghem (BCH) codes are utilized for SW coding. Also LDPC codes have been proposed and their performance has been analysed in [START_REF] Elkouss | Rate compatible protocol for information reconciliation: An application to QKD[END_REF]. However, as most of the work on coding for SW settings only consider long blocklengths, respective results for the short blocklength have not been reported.

In this paper we investigate coding for reconciliation where the two parties have generated a binary random sequence from their observations of a common random source. As an example, the possible source of random variables can be the wireless channel state or the challenge-response of physical unclonable functions (PUF). We implement and compare the reliability / complexity performance of polar, LDPC and BCH codes in the short blocklength regime. All the selected codes show promising performance in information transmission in the short blocklength.

The paper is organized as follows. In Section II, we describe the system model. In Section III, we explain the details of the codes' implementations, and in Section IV we present the numerical results and compare the performance of the codes investigated. Finally, we conclude the paper in Section V.

II. SYSTEM MODEL

In our system model, we assume Alice and Bob observe binary sequences Y A , Y B of length n by quantizing their local measurements of a common source of randomness, represented by X A and X B , respectively. At present, for the sake of simplicity, we assume that a passive adversary, referred to as Eve, is unable to obtain any information about the observed sequences.

We assume that the generated binary sequences are independent and identically distributed (i.i.d.) with bits at any index having equal probabilities of being 0 or 1, i.e.,

P r (Y

A [i] = 0) = P r (Y B [i] = 0) = 0.5. (1)
Also the corresponding bits of the sequences Y A , Y B are different with probability The code rate is defined as R = k n and the frame error rate (FER) is defined as the probability that Bob's estimation of Y A is erroneous P r(Y A = ŶA ). In this setup, Bob first estimates the sequence ŶA and then, Alice and Bob extract the common key K independently using privacy amplification. (Fig. 1). Motivated by the promising performance of LDPC [START_REF] Zhou | Improving short-length LDPC codes with a CRC and iterative ordered statistic decoding[END_REF], polar [START_REF] Tal | List decoding of polar codes[END_REF], and BCH codes [START_REF] Fossorier | Soft-decision decoding of linear block codes based on ordered statistics[END_REF] with list decoding in short blocklengths for information transformation applications, in this paper we implement and compare against the upper bound on SKG rates in [START_REF] Tyagi | Converses for secret key agreement and secure computing[END_REF], when these three families of SW decoders are employed.

p P r (Y A [i] = Y B [i]) = p, i = 1, . . . , n. (2)

III. CODES IMPLEMENTATION

A. LDPC codes with ordered statistic decoding LDPC codes are powerful error correcting codes that can approach the Shannon limit at very large blocklengths. They exhibit promising performance for a variety of communication channels [START_REF] Shakiba-Herfeh | LDPC code design for fast fading interference channels[END_REF]- [START_REF] Shakiba-Herfeh | LDPC codes for interference channels in the primary decodes cognitive regime[END_REF]. However, in general, they do not work well in short blocklengths. To address such shortcomings, LDPC codes enhanced with ordered statistic decoding (OSD) is one of the techniques proposed to achieve near maximum likelihood (ML) efficiency for short blocklengths [START_REF] Zhou | Improving short-length LDPC codes with a CRC and iterative ordered statistic decoding[END_REF]. The main concept behind OSD is that we start by selecting the k most reliable independent positions, where k is the code rank. Afterwards, we make hard decisions on the value of the selected bits, based on the log-likelihood ratios (LLR). Then, by flipping the values of up to t bits among them, we produce a candidate list of codewords. Finally, we choose the most likely codeword by performing an ML search in the list [START_REF] Wu | Soft-decision decoding using ordered recodings on the most reliable basis[END_REF]. The size of the candidate OSD list increases with respect to t as t i=0 k i . In our implementation, Alice sends her syndrome S = Y A H t , where H t is the transpose of the parity check matrix. At the receiver side, Bob first feeds Y B and S to the sumproduct (SP) algorithm to generate soft information of the bits (i.e., the LLRs). Then the LLRs are passed to the OSD block to estimate ŶA (Fig. 2). In the iterative version of LDPC code with OSD, to enlarge the candidate list of codewords, the LLRs after each iteration of SP algorithm are passed to OSD, and at the end the OSD block pick the most likely codeword. The size of the candidate list of codewords in iterative scheme is at most itr times larger than the non-iterative scheme, where the itr denotes the number of iterations of the SP algorithm.

B. Polar codes with list decoding

Polar codes are linear block error correcting codes that can provably achieve the capacity of a binary-input discrete memoryless channel as the code length goes to infinity. Using successive cancellation list decoding, which holds a list of most likely decoding paths, can significantly improve the performance of polar codes in finite blocklengths. The cyclic redundancy check (CRC) will boost the performance of list decoding even further. The CRC aids the decoder in selecting the correct decoding route from a list of options, even if it is not the most likely one. [START_REF] Tal | List decoding of polar codes[END_REF].

In our implementation (similar to [START_REF] Arikan | Source polarization[END_REF]), Alice encodes her sequence

Y A as U = Y A G n , where G n = 1 0 1 1
⊗n is the encoder matrix as defined in [START_REF] Arikan | Source polarization[END_REF]. Alice sends the syndrome S which contains S 1 and CRC bits with length l. S 1 has length n -k -l and contains high-entropy bits of U as follows

H(U [i]|Y A , U i-1 ) ≥ H(U [j]|Y A , U j-1 ), 1 ≤ i, j ≤ n, (3)
where i is the position of transmitted bits and H(•) denotes entropy. Therefore, the actual rate of the polar code is R = k+l n . On the other side, Bob applies CRC-aided successive cancellation list decoding to estimate ŶA . In the list decoding scheme, Bob tracks L decoding paths simultaneously, where L is the list size. Finally, the decoder picks the most likely codeword which also satisfies the CRC condition among the L paths (Fig. 3). Note that the complexity of list decoding polar coding grows linearly with the list size.

C. BCH codes with list decoding

BCH codes are a class of cyclic error-correcting codes constructed by polynomials over a finite field. One of the main features of BCH codes during code design is the number of guaranteed correctable error bits. A binary BCH code is defined by (n BCH , k BCH , t BCH ), where n BCH = 2 w -1 is the blocklength, k BCH is the message length, and t BCH represents the number of guaranteed correctable error bits. To improve their error correcting capability, BCH codes can be armed with list decoding.

In our implementation of list decoding, Alice calculates the syndromes as S = Y A H t , where H is the parity check matrix of the BCH code, and transmits it through the public channel. On the other side, Bob, first generates a candidate list by flipping up to t bits of the measured sequence Y B . After feeding the candidate list to the BCH decoder, the BCH decoder finds a possible solution for each element of the candidate list and picks the solution which is the most likely codeword with respect to the measured sequence Y B (Fig. 4). In this implementation of list decoding, the size of the list increases with respect to t as t i=0 n i . To reduce the size of list, without a significant degrade in the performance of decoder some heuristic methods exist in the literature [START_REF] Cho | Efficient software-based encoding and decoding of bch codes[END_REF].

IV. NUMERICAL RESULTS

In this section, we provide the FER performance of the aforementioned codes in the SKG setting for 128 and 512 bits. We also compare the FER performance of codes with the information theoretical finite length upper bound reported in [START_REF] Tyagi | Converses for secret key agreement and secure computing[END_REF], which translates to a lower bound in the FER.

For the first instance, we consider blocklength 128 bits. We implement a (n, k) = (128, 75) polar code with 11 bits CRC and a (127, 64, 10) BCH code. Also, we pick regular (3, 6) LDPC codes with length 128 bits 2 . The maximum number of iterations for SP algorithm is set to 50. In Fig. 5, the FER performances of half rate codes with 128 bits blocklength are depicted (i.e., the key length after privacy amplification is 64) and compared to the upper bound 2 We note that, the variable and check nodes degree distributions can be optimized to improve the performance of LDPC codes [START_REF] Shakiba-Herfeh | LDPC code design for fading interference channels[END_REF], however the degree optimization of the LDPC codes is out of the scope of this work.

reported in [START_REF] Tyagi | Converses for secret key agreement and secure computing[END_REF]. Note that the upper bound becomes inaccurate for short blocklength while their accuracy improves as the blocklength increases. That is the reason we observe comparatively a large gap between the upper bound and the codes' FER performances for n = 128 [START_REF] Hentilä | On polar coding for finite blocklength secret key generation over wireless channels[END_REF].

As it is depicted in Fig. 5, the polar code with list size 128, second order BCH list decoding code, and LDPC iterative list decoding for t = 2 outperform and show almost similar performance at FER 10 -3 . Moreover, The results show that, while traditional polar codes perform poorly in short blocklengths, their performance can be significantly improved by equipping them with list decoding. For example, the polar code with list size 128, provides around two orders lower FER compared to the classical polar code at H(Y A [i]|Y B [i]) = 0.1944. However, this improvement comes at the cost of 128 times more decoding complexity. We also observe that LDPC iterative list decoding outperforms the non-iterative one. Furthermore, list decoding also improves the performance of the BCH and LDPC codes, but the improvement with respect to the additional complexity is not as notable as in the case of the polar code.

Furthermore, to generate keys that can be used with standard block ciphers, e.g., AES-128 or AES-256, the SKG process is assumed to run with blocklength 512 and rate half. As it is demonstrated in Fig. 6, we observe a significant improvement in the performance of polar codes by using list decoding. On the other hand, we do not observe a noticeable improvement by using list decoding size t = 1 for LDPC and BCH code for this blocklength. We conclude this observations as the strength of list decoding with a fixed size as the blocklength increases, diminishes. We do not generate higher list decoding order for LDPC and BCH codes due to their high decoding complexity. Also, at blocklength n = 512, the gap between the FER of the polar code with list size 128 and the lower bound is less than the case n = 128. We posit that one of the reasons is that the lower bound at this length is tighter than the first instance.

In the simulations for short blocklength, we observe that although the considered implementation of list decoding for LDPC and BCH codes improves the performance of the code in short blocklength (∼100 bits). However, the improvement diminishes as we increase the blocklength, and to obtain enhancement in larger blocklength, we need to increase the order of the list decoding, which adds the complexity of decoders significantly. For these schemes utilizing heuristic methods, to modify the list size, can effectively reduce the complexity. On the other hand, we observe that polar codes with list decoding provide a significant improvement for medium size blocklength as well as short blocklength.

V. CONCLUSION

In this paper, we consider coding for Slepian-Wolf key reconciliation in the short blocklength regime. We assume that two legitimate parties observe binary sequences, generated from measurements of a common random source. We Fig. 6: Comparison of FER performance of codes with the lower bound in [START_REF] Tyagi | Converses for secret key agreement and secure computing[END_REF] for n = 512.

implement three different families of codes, including LDPC codes, BCH codes, and Polar codes to detect and correct the mismatches in the noisy observations. Our simulations for blocklengths of n = 128 and n = 512 bits, considering both reliability and decoding complexity, demonstrate that polar codes with a list size of 128 outperforms all other codes.

Fig. 1 :

 1 Fig. 1: The system model.

Fig. 2 :

 2 Fig. 2: Decoder structure for LDPC codes with ordered statistics decoding.

Fig. 3 :

 3 Fig. 3: Decoder structure for polar codes with list decoding.

Fig. 4 :

 4 Fig. 4: Decoder structure for BCH codes with list decoding.

  A [i] | Y B [i])=H(p) (bits)

  Fig.5: Comparison of FER performance of codes with the lower bound in[START_REF] Tyagi | Converses for secret key agreement and secure computing[END_REF] for n = 128.

	10	0
	10 -1
	FER									LDPC iterative list decoding t=2 LDPC list decoding t=2 LDPC list decoding t=1
										LPDC code
										LDPC iterative list decoding t=1
										Polar code
										Polar + CRC (list=8)
										Polar + CRC (list=32)
										Polar + CRC (list=128)
										BCH code
										BCH code list decoding t=1
										BCH code list decoding t=2
										Lower bound
	10 -1
	10 -2
	FER
	10	-3							LDPC code
										LDPC list decoding t=1
										Polar code
										Polar code + CRC (list=8)
	10 -4							Polar code + CRC (list=32) Polar code + CRC (list=128)
										BCH code
										BCH code list decoding t=1
										Lower bound
	10 -5	0.15	0.2	0.25	0.3	0.35	0.4	0.45	0.5
							H(Y A [i] | Y B [i])=H(p) (bits)

ACKNOWLEDGEMENTS

This work has been supported by the ANR PRCI project ELIOT (ANR-18-CE40-0030 / FAPESP 2018/12579-7) and the INEX Paris Universite projects PHEBE and eNiGMA.


Keqiang Li

Yajuan Cheng

Hongying Wang

Yangyu Guo

Zhongwei Zhang

Marc Bescond

Massahiro Nomura

Sebastian Volz

Xiaohong Zhang

Shiyun Xiong
email: xiongshiyun216@163.com

Phonon resonant effect in silicon membranes with different crystallographic orientations

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

INTRODUCTION

Designing structures with well controlled heat flow is extremely important in many fields ranging from waste energy harvesting, novel computations, information storage based on heat, to efficient heat spreading [1][2][3][4]. The phonon scattering method is a widely studied approach for heat transfer engineering, which either introduces scattering sources inside materials to impede phonon transport [2] or reduces structural disorder to facilitate phonon transport [5,6]. Although the phonon scattering method can tune the thermal conductivity (TC) of materials in a large magnitude, it relies on the phonon particle picture and is efficient in controlling high-frequency phonons but less efficient for lowfrequency long-wavelength modes. Moreover, the scattering sources are in general inside materials, which may affect other material properties, such as electrical conductivity. The second strategy to control phonon transport is the phononic crystal method [7,8], which is based on the fabrication of periodic holes inside materials. In phononic crystals, the interference between forward waves and back scattered waves modifies the band structure of phonons and create phonon band gaps. While similar to phonon scattering, the generated periodic holes inside materials could be also problematic to other properties. Recently, a new strategy for controlling phonon transport at nanoscale based on nanophononic metamaterial (NPM) has been proposed [9][10][11][12]. Similar to the phononic crystal approach, the NPM method is also based on the phonon wave effect by introducing phonon resonances. However, the resonant structures are usually located at the surface of a nanostructure, which is different from the phononic crystal method. As a result, the NPM strategy can be regarded as a complementary heat transfer engineering method to phonon scattering and phononic crystal methods. Isolated surface pillars were initially found to be ideal structures to create surface resonances [9,10,13,14]. It was found that the surface pillars can reduce the TC through the resonant hybridization effect, which can simultaneously decrease the group velocity and the relaxation time of phonons, especially at low frequencies [10,15,16]. The ability to tune the low-frequency phonon transport with small pillars makes the NPM an ideal approach to engineer the thermal transport at low temperatures, where only phonons with low energy are populated. Besides, it has been illustrated that phonon resonances can be complementary to phonon scattering for entire frequency phonon engineering, i.e., they can im-pede both low-and high-frequency phonon transport synergistically. Apart from pillared resonant structures, the phonon resonant effect has also been found in many other structures like nanojunctions [17,18], fullerene attached graphene [19], graphdiyne [20], clatherates [21][22][23][24]. Especially, the native oxide surface layers in Si-membranes were found to be able to generate random resonances [25], which can greatly simplify the NPM structure fabrications. Moreover, the machine learning strategy has been applied to optimize the distribution of resonant structures, which can help to further improve the thermal conductivity reduction magnitude by resonances. [26] More recently, the phonon resonant mechanism was further proved in graphene nanoribbon systems through an observation of unexpected TC enhancement by defects [27].

In this work, we aim to investigate the effect of phonon resonances on the thermal transport of Si membranes bounded with different crystalline facets based on nonequilibrium molecular dynamics simulations (NEMD). We uncover that the TC of Si membranes can be largely reduced by surface pillars. However, the reduction percentage strongly relies on membrane orientations, i.e., surface types. Besides, the quality of interfaces between surface pillars and the base material, which plays an important role in resonant hybridization strength, is further evaluated to control heat transport. We consider Si membranes bounded with {100}, {110} and {111} facets, for which the thicknesses are 4.3 nm, 3.8 nm and 3.8 nm, respectively; see Fig. 1(a). For the {100}-membrane, the two in-plane directions are equivalent ([100] direction) and simulations along one direction is considered. For the {110}-and {111}-membranes, the two in-plane directions are non-equivalent and the thermal transport along the [100] and [110] crystallographic directions for the {110}-membrane, and along the [11 2] and [1 10] crystallographic directions for the {111}membrane is studied. In cubic lattices, the lattice directions [11 2] and [1 10] are equivalent to [112] and [110], respectively. Hereafter we refer to [11 2] as [112] and to [1 10] as [110] for simplicity. Starting from the pristine membranes, we constructed the resonant structures by introducing pillars on the two surfaces; see Fig. 1(b). Two pillar heights, 2 nm and 8 nm, are considered to study the height effect on thermal transport. The periods of the resonant structures (Fig. 1b) are as follows: L 100 = 4.3 nm for the {100}-membrane, L 110 = 3.8 nm and L 100 = 4.3 nm for the {110}-membrane, L 112 = 4.0 nm and L 110 = 3.8 nm for the {111}-membrane. Here the subscript index denotes the crystalline orientation, which also corresponds to the transport direction. The side lengths of pillar cross-sections are one half of the structure period in the corresponding directions. To study the interfaces quality (between pillars and the base material) dependent phonon resonances, we further create vacancy defects (of size 1.1 × 1.1 × 1.1 nm 3 for {100}-membrane, 1.1×1.2×1.2 nm 3 for {110}-membrane and 0.7×1.2×0.9 nm 3 for {111}-membrane) around the interface between pillars and the membrane; see Fig 1(c). All TCs were calculated with NEMD simulations implemented in the GPUMD package [28] with the interatomic interactions described by the Stillinger-Weber three-body potential [29], which is commonly adopted for TC calculations in silicon-based materials [30][31][32]. A time step of 0.5 fs is chosen to ensure the energy conservation during the simulations. All simulations were performed at a target averaged temperature T = 300 K with the thermostat temperature T h/c = T ± ∆T /2, where ∆T = 40 K was used. The schematic illustration of NEMD set up is shown in Fig. 2. The atoms in the two ends (1 nm in length) are fixed to avoid the drift of samples. Close to the fixed regions, the atoms within the length L bath = 16 nm are coupled to the langevin thermostat [33,34] with a coupling constant of 0.25 ps. The width of samples (including the pristine structures) is set as one period length of the corresponding resonant structures. Simulations for various system lengths L sys are performed to extract the TC corresponding to the infinite length. In NEMD simulations, a temperature jump can usually be observed close to the thermostat due to the ballistic transport of a portion of phonons. It was proved by Li et al [35] that such a temperature jump in NEMD simulations is physical and should not be excluded for TC calculations. As a result, we calculated the TC based on the temperature difference ∆T according to the Fourier's law:

STRUCTURES AND METHODS

κ = J ∆T /L s , ( 1
)
where J is the heat flux and L s corresponds to the sample length. For both the pristine and resonant structures, the heat flux J is evaluated based on the cross-section of the base membrane. The obtained total heat current can be decomposed into each phonon frequency as proposed by Sääskilahti et al [36,37] and later modified by Fan et al [38]. The spectral heat current (SHC) is expressed as:

J(ω) = 2 ∞ -∞ dte iωt K T (t) (2)
where K T (t) is the component of K(t) in the transport direction. K(t) is evaluated as:

K(t) = i j =i r ij (0) ∂U j ∂r ji (0) • v i (t) (3)
here r ij is the distance vector between atoms i and j, U j denotes the potential energy of atom j and v i is the velocity vector of atom i. The summation is taken over all atoms in region R shown in Fig. 2. Here for samples with lengths shorter than ∼ 40 nm, the size of region R is slightly smaller than the system size to exclude the atoms connecting the thermal baths. For samples with lengths longer than ∼ 40 nm, the length of region R is set as ∼ 40 nm.

RESULTS AND DISCUSSIONS

Since the TC obtained directly from NEMD simulations corresponds to a finite sample length, one needs to extrapolate the TC to infinitely long structures. The linear extrapolation between 1/κ and 1/L as proposed by Schelling et al [39] has been widely applied to obtain the intrinsic TC due to its simplicity. The linear extrapolation was derived by assuming that all the phonon modes have the same averaged properties. However, it was recently pointed out that this method is only valid when the extrapolation is based on samples with lengths larger than phonon mean free path (MFP). Such a precondition dramatically increases simulation loads, especially for materials with long MFP. Actually, if the properties of each phonon mode are considered explicitly, the length dependent TC is a complex function of 1/L and can be expanded by a Taylor-series as shown by Sellan et al [40]. As a result, to improve the extrapolation accuracy, one can truncate the expansion to the second-order: [20,[40][41][42] 1

κ (L) = 1 κ 0 + α L + β L 2 (4)
where κ(L) is the TC obtained by NEMD simulations for samples with length L. κ 0 corresponds to the extrapolated TC with α and β the extrapolation coefficients. To improve the extrapolation accuracy, the values of short samples are omitted and the intrinsic TC is extrapolated with relatively long samples [40][41][42]. Fig. 3 reports the TC extrapolation for the pristine and resonant (2 nm in height) {110}-membrane based on Eq. 4. For both the pristine and resonant membranes, the relationship between 1/κ and 1/L clearly deviates from the linear shape, especially for pristine structures. In principle, when the majority of the phonons are travelling diffusively, i.e., for samples longer than phonon MFP, one could approximate the inverse of TC as a linear relationship with 1/L. However, for samples with lengths smaller than phonon MFP, i.e., the majority of the phonons are travelling ballistically, the obtained TC will deviate from the linear approximation. Such deviation could be compensated by adding a quadratic term of 1/L 2 in the TC extrapolation formula as shown in Eq. 4. In resonant membranes, the phonon MFP is dramatically reduced, especially at low frequencies (see below). As a result, their deviations from the linear relationship are much smaller than the case of pristine structures, where the phonon MFP is relatively large. It should be pointed out that, even for quadratic fitting, one should not use samples with ultra-short lengths [42]. In our case, the TC can be well extrapolated by the quadratic fitting for samples with lengths larger than 30 nm (Fig. 3). For pristine structures, since their phonon MFPs are relatively large, we performed simulations for long samples with maximum length reaching ∼ 550 nm to improve the extrapolation accuracy. Fig. 4(a) illustrates the TC along different crystallographic directions in the pristine and resonant membranes bounded with {100}, {110}, and {111} surfaces. The Stillinger-Weber potential adopted here overestimates the TC of silicon, however, we aim to study the phonon resonant effect on thermal transport in membranes with different orientations. Consequently, the relative TC change is more important and the absolute value of TC will not change our conclusions. Due to the symmetry breaking, the pristine {110}-membrane shows large anisotropic thermal transport, the TC along the [100] and [110] directions are 160 W/mK and 218 W/mK, respectively. Although the two in-plane transport directions are along different crystallographic directions in pristine {111}-membrane, it shows almost isotropic thermal transport with κ =∼ 225 W/mK. On the contrary, they are different in the case of the pristine 110membrane. This might be due to different bond symmetries for atoms on the {110} and {111} surfaces. For both the {110}-and {111}-membranes, the atoms on surfaces possess one dangling bond. However, among the three left bonds, two of them are on the surface and one is below the surface for the {110}-membrane. In contrast, the three left bonds in the {111}-membrane are all below the surface with equivalent environment (See Fig. 1(e)). As a result, a large/small anisotropic effect is expected for the {110}/{111}-membrane. The equivalent bond environment for surface atoms in the {111}-membrane can also lead to a small phonon-surface scattering rate, resulting in a larger TC compared to the {100}-and {110}membrane. With the introduction of resonant pillars on surfaces, the TC of all membranes is reduced greatly. However, the reduction percentage is dependent on the membrane type and transport direction as shown in Fig. 4(b). Among the three membrane types, the resonant effect is most effective in the {111}-membrane. In pristine membranes, the {111}-membrane possess the largest TC compared to the {100}-and {110}-membrane. However, with the creation of resonant pillars, the TCs along both the [110] and [112] directions are smaller than the values in the other two membrane types. Note that the pillar densities (the ratio between pillar and membrane length along a specific direction) along all directions in all membrane types are the same (50%). With resonant pillars of 2 nm in height, the reduction percentage can reach up to 81% and 83% along the [110] and [112] directions, respectively. The TC reduction in other membrane types is also larger than 70%, showing good TC tunability of phonon resonances. Interestingly, adding resonant pillars on the {111}-membrane will change the isotropic transport in pristine membranes to an anisotropic transport with the TC along the [112] direction smaller than the one along the [110] direction. Although the three bonding environment of surface atoms in the {111}-membrane is the same, the lattice constant is different along the [110] and [112] directions. As a result, the pillar cross-section lengths along the two directions are different, which eventually produces the anisotropic thermal transport. To evaluate the effect of pillar height on TC, the TC of resonant membranes with pillar height of 8 nm is also simulated. Pillars with increased height can lead to a further reduction of TC in all membrane types. However, such a further reduction is not significant compared to the original TC reduction by 2 nm height pillars.

To gain more knowledge on TC variations, we calculated the spectral heat current (SHC) in NEMD simula- S1. For pristine membranes, SHC at high frequencies decreases with the increase of sample length, which is a result of anharmonic scattering. However, at low frequencies, all structures show ballistic transport behavior with SHC almost unchanged with sample length up to ∼ 500 nm. The frequency range of ballistic transport depends on the system. To compare quantitatively, we define a critical frequency ω c , below which SHC decreases less than 10% when the sample length is increased from the shortest one to the longest one as demonstrated in Fig. 5(a)-(c). The corresponding frequency ranges are demonstrated by the pink regions in Fig. 5(a)-(c). We find that the critical frequency depends strongly on transport orientation while less affected by surface types. The obtained critical frequencies for the {100}-[100], {110}-[100], {110}-[110], {111}-[110], and {111}-[112] structures are 1.4, 1.5, 2.0, 2.0, and 1.9 THz, respectively. The ballistic transport region is smaller in the transport direction [100] compared to the [110] and [112] directions, no matter the surfaces are bounded by {100} or {110} facets. The relatively smaller ballistic transport region can result in smaller TC along [100] direction. Although the surface orienta-tion does not have a large impact on the frequency range of ballistic transport, it can affect SHC at higher frequencies. For example, the SHC in {110}-[100] beyond ω c is smaller than that in {100}-[100] with similar sample lengths (Fig. S1). Together with the small ω c , the {110}-[100] structure eventually results in the smallest TC among the five studied cases. In the presence of surface pillars, the SHC is reduced noticeably compared to the corresponding pristine structure of similar length. Especially at long sample lengths, SHC is dramatically reduced. More importantly, SHC is notably reduced with the increase of sample length at low frequencies, where ballistic transport appeared in the corresponding pristine structures. Such a phenomenon indicates that the surface resonant pillars can dramatically tune the low-frequency thermal transport in all membranes, which is in good agreement with previous findings [10,25,32]. The ability of tuning the low-frequency phonon transport is an advantage of the resonant mechanism over scattering, which usually impedes high-frequency modes extensively [10,25]. As observed from TC calculations, although the TC in the pristine {110}-[110] and {111}-[112] structures is similar, adding surface pillars can make their TC much different. From SHC curves (Fig. 5(e)-(f)), one could observe that the small TC in the resonant {111}-[112] structure is originated from the frequency range between 3 and 12 THz, where the SHC is smaller in the {111}-[112] membrane compared to the {110}-[110] membrane of similar length. The TC along the [100] direction of resonant {110}-membrane is also small, however, different from the case {111}-[112], it possesses small SHC beyond 12 THz (Fig. S1).

The resonant effect relies on the hybridization between resonant and propagating modes. As a result, the interface quality between the base membrane and pillars is crucially important as it determines the penetration of resonant phonons into the base membrane. To illustrate such an effect, we examined the TC of {100}-membrane with pillars in [110] and [111] crystallographic orientation of 2 nm in height. In such cases, the interface between the {100}-membrane and [110]/[111] pillars is not perfect anymore. The obtained TCs for pillars in [110] and [111] orientation are 54.0±1.6 and 60.1±2.6 W/mK, respectively. These values are larger than the membrane with [100] pillars (51.3±2.1 W/mK). The larger TC in 100-membrane with non-perfect interfaces ([110] and [111] pillars) reveals that the non-perfect interface can reduce the resonant hybridization. Actually, it was shown by Wang et al that the TC of resonant structure can be even enhanced by introducing vacancies at the interface between resonant pillars and the base material compared to the perfect resonant structures [27]. Such an anomalous TC enhancement is also due to the decoupling between resonant and propagating phonons. To further investigate such decoupling effects by interface qualities in our membranes with different orientations, we also created vacancies inside the material as demonstrated in Section II. Fig. 6 compares the TC of resonant structures with and without interface vacancies. For all structures, an increase of TC is observed with the introduction of interface vacancies when compared to the perfect resonant structures, which inhibits the penetration of resonant modes into the base material. The TC enhancement is most noticeable in the {111}-membrane, where the enhancement along the [110] and [112] directions reaches 18% and 20%, respectively. The large enhancement of TC in the {111}-membrane is in accordance with the significant resonant effect in this system. In contrary, the TC enhancement by interface vacancies is small in the {100}-membrane and along the [110] direction of the {110}-membrane due to the relatively weak resonant effect. However, we found that the TC enhancement effect is also large along the [100] direction in the {110}-membrane, which possesses the smallest resonant effect. Such an uncorrelated TC enhancement with resonant strength indicates that the resonant coupling in the {110}-[100] structure is very sensitive to the interface quality. The enhancement of TC by interface vacancies not only proves the existence of phonon resonant hybridization in pillared structures, but also brings a degree of freedom to tune the TC in resonant structures.

To gain more knowledge on the frequency range where heat flux is primarily affected by the interface vacancies, the SHC is analyzed for the longest samples that we studied. Three systems, {100}-[100], {110}-[110], and {111}-[112], are selected for analysis and the corresponding results are illustrated in Fig. 7. The SHC of the other two cases, {110}-[100] and {111}-[110], is provided in Fig. S2. Interestingly, for all structures, the SHC is only enhanced in a narrow frequency window, not at all frequencies. The range of frequency window in which SHC is enhanced depends on transport directions. As illustrated by Fig. 7(a) and Fig. S2(a), SHC is enhanced between 5 and 9 THz when the transport direction is along the [100] crystallographic direction, no matter it is in the {100}-or {110}-membrane. This frequency range is situated immediately after the first peak of the SHC. When the heat flux is along the [110] crystallographic direction, SHC is enhanced around the first peak in both the {110}-and {111}-membranes (Fig. S2(b)). The frequency range is also enlarged to [4][5][6][7][8][9] THz. For {111}-[112], SHC is also enhanced around the first peak but the frequency range is further broadened to [4][5][6][7][8][9][10][11] THz compared to the transport along the [110] direction. As a result, one can conclude that the interface vacancy in resonant structures mostly decouples the hybridization between resonant modes and propagating modes at intermediate frequencies.

CONCLUSIONS

Based on NEMD simulations, the thermal transport along different crystallographic orientations in silicon membranes bounded with {100}, {110} and {111} facets is evaluated. Surface pillars are introduced to create resonant modes, which can tune the thermal transport in nanostructures with a large degree of freedom. It was found that in pristine membranes, the crystal symmetry breaking of bulk silicon by surfaces leads to anisotropic transport. SHC analysis demonstrates that there is a ballistic transport region at low frequencies in the pristine membranes. The frequency range of ballistic transport is mainly dependent on the transport direction while less sensitive to the surface orientation. The introduction of resonant pillars can effectively break the ballistic transport characteristics at low-frequencies, indicating the ability of tuning the low-frequency phonon transport. The reduction of TC by resonances strongly depends on both the transport direction and surface orientation. The resonant effect is maximized in the {111}membranes, where the TC is the largest and smallest before and after the introduction of resonant pillars among the three types of membranes. Along both the [110] and [112] directions of the {111}-membrane, the TC reduction by resonant pillars (2nm) reaches up to 83%. Along the [100] direction, the TC reduction is the smallest with the reduction percentage around 74% in both the {100}and {110}-membranes. Introducing vacancies at the interfaces between resonant pillars and the base material can reduce the resonant phonon hybridization by affecting heat transport at intermediate frequency range. The ability of tuning the resonant coupling by interface quality thus provides a further degree of freedom to engineer the resonant effect and thermal transport.

FIG. 1 .

 1 FIG. 1. Schematic illustration of (a) pristine Si membrane, (b) resonant Si membrane, (c) resonant Si membrane with interface vacancies, (d) cross-section of pillars indicated by the dashed line in (c) (top) and the top view of resonant structures (bottom), and (e) relative orientations between the four bonds of a Silicon atom and the three facets. The coordinates below indicates the crystallographic orientations of the {100}-membrane, {110}-membrane and {111}-membrane, respectively. All membranes lie in the x-y plane.

FIG. 2 .

 2 FIG.2. The set-up of NEMD simulations. Region R is defined for SHC calculations.

FIG. 3 .

 3 FIG.3. The extrapolation of TC based on Eq. 4 for pristine and resonant (2 nm) {110}-membrane. The size of error bars for all data is smaller than the size of the symbols.

FIG. 4 .

 4 FIG. 4. (a) The TC along different directions in pristine and resonant membranes bounded with different facets. (b) The TC reduction percentage by resonant pillars with height of 2 nm and 8 nm.

FIG. 5 .

 5 FIG. 5. The spectral heat current in NEMD simulations for pristine and resonant (2 nm) membranes with different sample lengths. (a) pristine {100}-[100], (b) pristine {110}-[110], (c) pristine {111}-[112], (d) resonant {100}-[100], (e) resonant {110}-[110], and (f) resonant {111}-[112]. For pristine membranes, the pink region denote the ballistic transport region with SHC reduces less than 10% from the shortest sample to the longest sample calculated.

FIG. 6 .

 6 FIG.6. Comparison of TC between perfect resonant structures and resonant structures with interface vacancy.

FIG. 7 .

 7 FIG. 7. Comparisons of spectral heat current before and after the introduction of interface vacancies for (a) {100}-[100], (b) {110}-[110], and (c) {111}-[112] structures. The length of the samples is indicated in the corresponding figures.

ACKNOWLEDGEMENTS

This work was supported by the National Natural Science Foundation of China under Grant No. 11804242, the CREST JST under Grant No. JPMJCR19I1, the Major Research Plan of the National Natural Science Foundation of China (Grant No. 91833303), and the Major International (Regional) Joint Research Project of the National Natural Science Foundation of China (Grant No. 51920105005).


Diane Coursier

David Coulette

Hélène Leman

Emmanuel Grenier

Gabriel Ichim
email: gabriel.ichim@lyon.unicancer.fr

Gabriel Ichim Live

Live-cell imaging and mathematical analysis of the "community effect" in apoptosis

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

Introduction

Inconspicuous at the single cell level, apoptosis is a powerful mechanism when triggered at the organismal level. It is the sculpting force behind all multicellular organism morphogenesis, ensuring all organs have the right size, supernumerary cells are removed and neurons innervate the correct targets.

In cancer, apoptosis is an efficient guardian against oncogenic transformation and it is the preferred mechanism employed by tumor suppressor proteins to eliminate dangerous rogue cells [START_REF] Strasser | Cell Death in the Origin and Treatment of Cancer[END_REF]. For instance, p53 can trigger apoptosis by up-regulating the expression of several pro-apoptotic proteins [START_REF] Strasser | Cell Death in the Origin and Treatment of Cancer[END_REF]. Apoptosis is undeniably the main executioner of most cancer drugs and radiotherapy [START_REF] Carneiro | Targeting apoptosis in cancer therapy[END_REF].

Apoptosis is triggered when one of two distinct pathways is engaged: the extrinsic or the intrinsic. Extrinsic apoptosis relies on the stimulation of a death receptor-family member (e.g., TRIAL-R, FAS or TNF-R) by a TNFrelated cytokine, and activation of caspase 8 and 10 via the death-inducing signaling complex (DISC) [START_REF] Dickens | The 'complexities' of life and death: death receptor signalling platforms[END_REF]. The intrinsic pathway is initiated by an intracellular death stimulus such as DNA damage, chemo-and radiotherapy or tumor suppressor activation. In this case, the point of no return is mitochondrial outer-membrane permeabilization (or MOMP) followed by cytochrome c release, apoptosome assembly and activation of caspase 9 and then 3 and 7. Mitochondrial permeabilization is tightly regulated by the BCL-2 family proteins that can be anti-apoptotic (BCL-2, BCL-xL or MCL-1), proapoptotic (BID, BIM, BAD, PUMA or NOXA) and effector pore-forming proteins such as BAX or BAK (for review see [START_REF] Bock | Mitochondria as multifaceted regulators of cell death[END_REF]).

Once apoptosis is triggered, this does not inevitably lead to cell death, as it is now accepted that cells can survive if mitochondrial permeabilization does not affect all mitochondria or if the ESCRT-III complex restores plasma membrane integrity (5) (6) [START_REF] Gong | To the edge of cell death and back[END_REF]. Currently, there is a lack of understanding of how and why the kinetics of apoptosis varies so much between cells or following different treatments.

The heterogeneous response to apoptosis is also described in isogenic cancer cell lines. This is puzzling since these cells are genetically identical and should express similar levels of both pro-and anti-apoptotic proteins.

Cellular heterogeneity to undergo apoptosis thus appears to originate from non-genetic transcriptional variability. This translates into cell-to-cell variations of pro-and anti-apoptotic protein expression and activity. For instance, the Lahav group elegantly showed that there is a high variability of p53 response to DNA damaging agents; apoptotic cells accumulate p53 much faster and earlier, while this accumulation is delayed in resistant cells, allowing time for the up-regulation of pro-survival IAP (inhibitor of apoptosis) proteins [START_REF] Paek | Cell-to-Cell Variation in p53 Dynamics Leads to Fractional Killing[END_REF].

Recent publications shed light on the effects apoptotic cells exert on neighboring cells. They can actively release EGFR ligands, FGF2 or Wnt3 to boost the survival and proliferation of neighboring cells (9) [START_REF] Bock | Apoptosis-induced FGF signalling promotes non-cell autonomous resistance to cell death[END_REF] [START_REF] Ankawa | Apoptotic cells represent a dynamic stem cell niche governing proliferation and tissue regeneration[END_REF]. This paracrine crosstalk is also complemented by mechanotransduction mechanisms triggered by tissue stretching around apoptotic cells, which involves the Yes-associated protein (YAP) pathway [START_REF] Kawaue | Mechanics defines the spatial pattern of compensatory proliferation[END_REF].

There is therefore an increasing need to accurately estimate the variability of apoptosis induction in a given cell population. To achieve this, we developed a medium-throughput imaging pipeline using apoptotic cell markers, coupled with mathematical analysis. This allowed us to unveil a neighboring effect for apoptosis induction, irrespective of cell death stimuli or cell density, whereby the probability of undergoing apoptosis increases in the proximity of apoptotic cells. Hence, the combined application of imaging and computational analysis to evaluate the response to a given apoptotic stimulus may provide new insight into misunderstood phenomena such as fractional killing or apoptosis-induced proliferation, which are major concerns in the context of apoptosis-based cancer therapies.

Results

Live-cell imaging system to trigger and measure multi-stimuli apoptosis

To mathematically analyze the kinetics of apoptotic cell death while varying different parameters, we first set up a cellular system allowing medium-throughput quantification of apoptosis at the cell population level. For this, we used melanoma WM115 cells stably expressing mCherry-tagged histone H2B (H2B-mCherry) to spatially identify the exact topology of a given cell (Figure 1A). Following treatment with an apoptotic stimulus, the induction of apoptosis is easily visualized using the live cell-impermeant dye SYTOX Green (SG). Live cells can thus be identified by the nuclear-localized H2B-mCherry-only signal (red nuclei), whereas apoptotic cells appear yellow, owing to the colocalization of H2B-mCherry and green SG (Figure 1A). Next, we used various stimuli to trigger apoptosis. Treatment with TNFα and cycloheximide (CHX) was initially used as a model for death receptormediated apoptosis. Of note, CHX enhances TNFα-induced apoptosis by blocking the translation of short-lived anti-apoptotic proteins [START_REF] Ruff | Rabbit tumor necrosis factor: mechanism of action[END_REF]. Downstream of the TNFα receptor 1/2 (TNFR1/2), the DISC activates caspase 8 and triggers MOMP via BID cleavage (Figure 1B). To induce intrinsic apoptosis, we used two death stimuli: first, a combination of BH3 mimetics ABT-737 inactivating the anti-apoptotic proteins BCL2, BCL-xL and BCL-w, and UMI-77 targeting MCL1; second, a doxycycline (dox)-inducible WM115 cell line to activate the mitochondrial pore-forming protein BAX (Figure 1B). These apoptotic stimuli were then applied to SG-pre-loaded WM115 H2B-mCherry cells plated at different densities and imaged periodically (1 or 2 hours intervals) using the IncuCyte live-cell imager (Figure 1A).

Cell clustering has no effect on the incidence of apoptosis

Following different pro-apoptotic treatments (Figure 1B) and IncuCytebased live cell imaging, we obtained datasets of H2B-mCherry images (red images) and the corresponding SYTOX Green staining (green images) at different time points, for a kinetics analysis of 48 hours. We then designed two independent pipelines to segment H2B-mCherry-marked nuclei, one based on the CellProfiler software and the other written in Python, using standard image libraries to ensure a robust analysis (Figure 1A). The analysis pipeline included the following steps. First, for a pair of red/green images the output is a segmented "red" image with a known red object centroid position (X i , Y i ) onto which a measurement of the green signal is superposed (I i , representing the maxima green signal intensity) (Figure 3A). Second, for a given time-point t, the analysis pipeline retrieves a set of N red centroid positions (X i , Y i ) with their associated SYTOX Green intensity signal I i . Data analysis is then performed using three parameters δ, n and I t , where δ is the inter-cellular distance at which another cell(s) is/are considered to be in contact and could influence each other, n is the threshold number of neighboring cells above which a cell is considered to be in a cluster and I t is the threshold intensity of SYTOX Green, derived from background noise, above which the cell is considered to be apoptotic. Of note, two cells are considered neighbors if the distance between their centroids is bellow δ (Figure 3A).

Once we determined δ, n and I t , we determined for each cell whether it was apoptotic (A) or not (Ᾱ). We then established the relationships of proximity between cells (the so-called neighboring effect) by determining for each cell whether it was in a cluster of cells regardless of their apoptotic status (C) or with a known apoptotic status (Ca) (Figure 3B). This method is well suited to quantify spatial relations between apoptotic cells.

Using this approach, we were first interested in determining whether a correlation between cell clustering and apoptosis (between C and A). This is driven by previous work suggesting that cell density is an important modulator of both extrinsic and intrinsic apoptosis [START_REF] Reuven | The Hippo pathway kinase Lats2 prevents DNA damage-induced apoptosis through inhibition of the tyrosine kinase c-Abl[END_REF] [START_REF] Bar | Attenuation of the p53 response to DNA damage by high cell density[END_REF][START_REF] Ma | Cell density plays a critical role in ex vivo expansion of T cells for adoptive immunotherapy[END_REF][START_REF] Kühn | Cell density modulates apoptosis in human articular chondrocytes[END_REF].

By examining three scenarios, in which a given cell was considered to be in a cluster if it was surrounded by 1, 2 or 3 neighboring cells (regardless of their apoptotic status), we found that a cell with no neighbors (C) had the same probability of entering apoptosis as a cell in a cluster (of 1, 2, or 3 neighboring cells) (Figure 3C-E for TNFα/CHX treatment). We then computed and plotted the time dependent p-value of a Fisher exact test that evaluates the independence between A and C (Figure 3C-E, bottom panels). Of note, if the p-value is high (as in Figure 3C-E, bottom panels), this means that the hypothesis that A and C are independent cannot be rejected. In other words, this means that the probability of being apoptotic whenever one cell has neighboring relation(s) (probability p(A/C), blue curve in Fig. 3C-E, top panels) and the probability of being apoptotic when one cell is isolated (probability p(A/C), orange curve) might be equal.

The same holds true when triggering intrinsic apoptosis with BH3 mimetics (Figure 3F) or when overexpressing the pro-apoptotic protein BAX (Figure 3G). Though these results were obtained at low cell density, we confirmed their relevance for cells grown at a high density (Supplementary Figure 1).

Hence, these data suggest that cell clustering has no impact on the occurrence of apoptosis.

The induction of apoptosis is enhanced by the proximity of a cell to an apoptotic cell

Next, we tested whether the proximity of a cell with one apoptotic cell or an apoptotic cluster (up to three apoptotic neighboring cells) was correlated with its propensity to undergo apoptosis. We obtained results in complete contrast to Figure 3, as the probability to enter apoptosis when in proximity with one or more apoptotic cells (probability p(A/C a ), blue curve) was significantly higher than the probability to enter apoptosis when isolated from apoptotic cells (probability p(A/C a ), orange curve) (Figure 4A-C for TNFα/CHX treatment). Indeed, the p-value (Fisher's exact test) was significantly low (p < 0.05) for the intermediate time points, meaning that both probabilities cannot be considered as equal. We observed a higher p-value at the beginning and the end of the kinetics, which does not affect our analysis, since there are few apoptotic cells at beginning of the time-course and, at the end, most cells are apoptotic and the neighboring effect saturates. This observation also holds true for treatment with BH3 mimetics and BAX overexpression (Figure 4 D,E). Finally, we validated these data at a higher cellular density (Supplementary Figure 2).

Taken together, our analysis pipeline reveals that apoptosis induction is influenced by a neighboring effect: the likelihood of a cell to undergo apoptosis is significantly higher when it is in proximity with apoptotic cells.

Discussion

Defective induction of cell death is the cause on many diseases, ranging from cancer to neurodegenerative diseases, and it is therefore important to study its propagative nature using appropriate tools.

Here, we report the development of an analysis pipeline for mediumthroughput imaging of apoptotic cancer cells that integrates several key parameters: temporal dimension, topographic localization of both living (cells stably expressing the mCherry-tagged histone 2B) and apoptotic (SYTOX Green-positive) cells, while applying several apoptotic stimuli.

Using this approach, we found that cancer cells have a similar probability of undergoing apoptosis if they are in proximity or not with other cells of unknown status. In other words, cell crowding does not have a protective effect against apoptosis. However, a cell is more likely to become apoptotic if it is in close contact with apoptotic neighboring cells. This holds true for all of the apoptotic stimuli we tested, either extrinsic or intrinsic inducers of apoptosis. These results are in line with Bhola and Simon's study demonstrating that daughter cells, with a high likelihood of being in proximity, synchronously undergo apoptosis. This synchrony is lost as cells become less related over time [START_REF] Bhola | Determinism and divergence of apoptosis susceptibility in mammalian cells[END_REF]. Apoptotic synchrony is therefore transiently heritable and is lost due to noise in protein translation [START_REF] Spencer | Nongenetic origins of cell-to-cell variability in TRAIL-induced apoptosis[END_REF]. To take into account progeny and monitor cell proliferation history, we would need to introduce into our imaging pipeline a lineage tracer such as the one developed by Oren and colleagues [START_REF] Oren | Cycling cancer persister cells arise from lineages with distinct programs[END_REF].

Our experimental settings are in 2D cellular cultures, yet one might apply it to 3D cultures such as cancer organoids. For this, however, the sample will require transparization and 3D imaging microscopy. Here, the proof-of-concept is provided for apoptosis, though the protocol could easily be adapted to other types of programmed cell death, such as caspaseindependent cell death, necroptosis or pyroptosis. Importantly, this method could also be helpful in identifying factors facilitating or inhibiting cell death propagation. Similar methods were developed to quantify the occurrence of fractional killing, a barrier for effective cancer treatments, and helped establish mediators such as the anti-apoptotic protein MCL1 (21) [START_REF] Inde | Kinetic Heterogeneity of Cancer Cell Fractional Killing[END_REF].

Recent research clearly demonstrates that while apoptosis is immunologically silent, it still involves releasing a plethora of signaling molecules that impinge on the fate of surrounding cells. The term of apoptosis-induced apoptosis was coined in Drosophila, where apoptotic cells in the wing disk triggered apoptosis in trans by secreting Eiger (the TNF homologue in flies), and stimulating pro-apoptotic JNK signaling [START_REF] Pérez-Garijo | Apoptotic cells can induce nonautonomous apoptosis through the TNF pathway[END_REF]. The same holds true in mammals during the regressive phase of the hair follicle [START_REF] Pérez-Garijo | Apoptotic cells can induce nonautonomous apoptosis through the TNF pathway[END_REF]. In cancer, the radiation-induced bystander effect, through which irradiated cells lead to the death of non-irradiated neighboring cells, is also a good example of this concept [START_REF] Riegman | Population Dynamics in Cell Death: Mechanisms of Propagation[END_REF]. Interestingly, recent studies revealed the rapid propagation of ferroptotic cell death, and the analysis pipeline described herein might provide further mechanistic insights [START_REF] Linkermann | Synchronized renal tubular cell death involves ferroptosis[END_REF] [START_REF] Riegman | Ferroptosis occurs through an osmotic mechanism and propagates independently of cell rupture[END_REF]. Inversely, apoptosis-induced proliferation has been unveiled and describes a state in which dying cells release prostaglandins and Wnt3, instructing the proliferation of nearby cells [START_REF] Li | Apoptotic cells activate the "phoenix rising" pathway to promote wound healing and tissue regeneration[END_REF] [START_REF] Ankawa | Apoptotic cells represent a dynamic stem cell niche governing proliferation and tissue regeneration[END_REF]. Recently, the group of Tait showed that apoptotic cells also release FGF2 to enhance the anti-apoptotic resilience in trans [START_REF] Bock | Apoptosis-induced FGF signalling promotes non-cell autonomous resistance to cell death[END_REF].

The analysis method developed here could increase and enrich our knowledge on "community effects" of cell death so that they can be blocked or enhanced, according to the physio-pathological context.

Material and methods

Cell lines

Human melanoma cells WM115 (a gift from R. Insall -The Beatson Institute, Glasgow, UK) were maintained in DMEM supplemented with 2 mM Lglutamine (ThermoFisher Scientific, 25030-24), non-essential amino acids (ThermoFisher Scientific, 11140-035), 1 mM sodium pyruvate (ThermoFisher Scientific, 11360-039), 10% FBS (Eurobio, CVFSVF00-01) and 1% penicillin/streptomycin (ThermoFisher Scientific, 15140-122).

Stable cell line generation

For retroviral transduction, Phoenix Ampho 293T cells (2x10 6 in 10 cm in diameter Petri dishes) were transfected with pQXIN H2B-mCherry (gift from Dr. Luca Fava) using Lipofectamine 2000 (ThermoFisher Scientific, 11668019). Twenty-four hours later, the virus-containing supernatant was filtered and used to infect WM115 cells in the presence of 1 µg/mL polybrene (Sigma-Aldrich, H9268). This was repeated with fresh viral supernatant 48 hours later. Melanoma cells stably expressing H2B-mCherry were further purified by cell sorting. WM115 H2B mCherry cells with doxycycline-inducible expression of the BAX protein were created using the Sleeping Beauty transposon system [START_REF] Kowarz | Optimized Sleeping Beauty transposons rapidly generate stable transgenic cell lines[END_REF]. The plasmids used were piTR1 BAX (gift from David Goldschneider) and pCMV(CAT)T7-SB100 (Addgene, 34879) for the expression of SB100X transposase. WM115 H2B mCherry cells were then cotransfect with piTR1 BAX and pCMV(CAT)T7-SB100 plasmids using Lipofectamine 2000 and selected with puromycin (1 µg/mL).

Western blotting

Cell lysates were prepared using NP-40 lysis buffer (1% NP-40, The nitrocellulose membranes were rinsed 3 times for 10 min in TBS-Tween 0.1% and then incubated with appropriate secondary antibody coupled to HRP for 1 hour at room temperature. Finally, the proteins of interest were detected using Clarity Western ECL (Biorad, 1705060) and chemiDoc imager (Biorad, 17001401).

Incucyte imager-based cell viability assay

Cell viability was determined using an IncuCyte Zoom imaging system (Sartorius). Cells were seeded onto Imagelock 96-well plates in media containing 30 nM SYTOX Green (Life Technologies, S7020). Following different apoptotic treatments, the cells were imaged every 60 or 120 minutes and the analysis of the SYTOX Green-positive cells was performed using the available IncuCyte image analysis software (Essen Bioscience).

Modeling the progression of apoptosis

Data description

Nuclei, set as objects, are delineated at each time step using automatic segmentation of mCherry fluorescence ("Red") images via the CellProfiler pipeline for primary objects, which is tuned to detect blob-like structures. The "apoptosis" signal is obtained by measuring the maximum intensity over each nucleus object of the "Green" fluorescence image. A threshold value for this intensity I t is chosen, which defines the fluorescence intensity considered as sufficient for a given nucleus to be considered as apoptotic. In the results presented herein, the threshold chosen was 0.2 that corresponds to 2/3 to 1/2 of the maximum intensity of the Green signal.

Statistical test

For a given time frame, segmentation data consist in a sequence of nuclei positions (X i ,Y i ) and green signal intensity I i , with i varying from 1 to N, where N is the total number of nuclei objects for this time frame.

Apoptotic variable A

Given the intensity threshold I t , we defined the class A of apoptotic cells as the set of nuclei such that I i ≥ I t .

Clustering (density) C and apoptotic clustering Ca variables

In order to separate "clustered" cells from isolated ones, we set a characteristic interaction distance δ and a threshold n on the number of neighbors. An object at position (X i ,Y i ) is considered to be part of a cluster (i.e., in the class C) if there are a least n nuclei inside the disk of radius δ and center (X i ,Y i ). An object at position (X i ,Y i ) is considered to be part of an apoptotic cluster (i.e., in the class C a ) if there are a least n apoptotic nuclei (i.e., in the class A) inside the disk of radius δ and center (X i ,Y i ).

Independence tests

We used classical Fisher's exact test to determine if random variables A and C are independent and subsequently if random variables A and C a are independent. In this setting, for the simplicity of notation, we identify the class A with the random variable that indicates if an object belongs or not to the class A, and similarly for classes C and C a . We gave exact p-values for this test. Generally, the null hypothesis (the considered pair of random variables are independent of each other) can be rejected when the p-value is lower than 0.05.

  Next, we validated the induction of apoptosis by assessing the processing of effector caspase-3 into p17/p19 fragments and the cleavage of PARP1, which is a hallmark of apoptosis effectiveness. As shown In Figure2A and B, treatment with TNFα/CHX (TC), BH3 mimetics (BH3m) and increasing doses of doxycycline efficiently induced apoptosis. Of note, caspase-3 processing and PARP1 cleavage were blocked by treatment with the pan-caspase inhibitor Q-VD-OPh. In addition, IncuCyte-based live-cell imaging was used to delineate the kinetics of apoptosis triggered by the different stimuli described above (Figure 2C, D). Apoptotic cells were marked by SG and quantified over time, with image acquisitions every 1-2 hours over a period of 48 hours on average. Figure 2E shows this imaging method can efficiently detect and discriminate between living cells (mCherry-positive) and apoptotic cells (SG-positive). The figure presents results of TNFα/CHX and BH3 mimetics treatment, albeit the same was observed for WM115 cells overexpressing BAX (data not shown). Using this cell model of apoptosis induction, coupled with mediumthroughput live-cell imaging we therefore established an imaging technique able to monitor the heterogeneous induction of apoptosis at the cell population level.

Figure 1 .

 1 Figure 1. Imaging and quantification of apoptosis induction using live and dead cell markers. A. Overview of the live/dead cell imaging and computer analysis protocol. B. Apoptotic signaling induced by treatment with TNFα/CHX or BH3 mimetics, or through BAX overexpression (aCasp8 -active Caspase-8; aCasp 3/7active Caspase-3 and -7, tBID -truncated BID).

Figure 2 .

 2 Figure 2. Validation of the apoptosis induction protocols. A. WM115 H2B-mCherry cells were treated either with TNFα/CHX (TC) (50 ng/mL of TNFα and 5 µg/mL CHX) or BH3 mimetics ABT-737/UMI-77 (10 µM each) for 12 hours in the presence or absence of the pan-caspase inhibitor Q-

Figure 3 .

 3 Figure 3. Apoptosis occurs independently of cell clustering. A-B. Definition of the parameters taken into account for the mathematical analysis of the effect cell clustering has on apoptosis induction. Cells with the nucleus half in red half in green are either alive or apoptotic, respectively. C-E. Probability estimates (upper panels) and the corresponding p-value (lower panels) for the occurrence of apoptosis in TNFα/CHX-treated WM115 H2B-mCherry cells surrounded by 1, 2 or 3 neighboring cells. F-G. Similar analysis as in (C-E), for WM115 H2B-mCherry cells treated with ABT-737/UMI-77 (F) and doxycycline (G). The neighboring effect for 1 cell is shown.

Figure 4 .

 4 Figure 4. Apoptotic cell clusters influence the life-death decisions in neighboring cells. A-C. Probability estimates (upper panels) and the corresponding p-value (lower panels) for the occurrence of apoptosis in TNFα/CHX-treated WM115 H2B-mCherry cells surrounded by 1, 2 or 3 apoptotic neighboring cells. D-E. Similar analysis as in (A-C), for WM115 H2B-mCherry cells treated with ABT-737/UMI-77 (D) and doxycycline (E). The neighboring effect for 1 apoptotic cell is shown.

Acknowledgements

Funding from Institute Convergence PLAsCAN (ANR-17-CONV-0002), LabEx DEVweCAN (University of Lyon), Agence Nationale de la Recherche (ANR) Young Researchers Project (ANR-18-CE13-0005-01), La Ligue Nationale Contre le Cancer and Fondation de France supported this work. We thank Brigitte Manship for reviewing the manuscript.

Author contributions


Gustavo A Nunez Segura

Arsenia Chorti
email: arsenia.chorti@ensea.fr

Cintia Borges
email: cintia@usp.br

Distributed DoS Attack Detection in SDN: Tradeoffs in Resource Constrained Wireless Networks

Keywords: Software-defined networking, intrusion detection, wireless sensor networks

The Software-defined networking (SDN) paradigm centralizes control decisions to improve programmability and simplify network management. However, this centralization turns the network vulnerable to denial of service (DoS) attacks, and in the case of resource constrained networks, the vulnerabilities escalate. The main shortcoming in current security solutions is the tradeoff between detection rate and complexity. In this work, we propose a DoS attack detection algorithm for SDN resource constrained networks, based on recent results on non-parametric real-time change point detection, and lightweight enough to run on individual resource constrained devices. Our experiment results show detection rates and attacker identification probabilities equal or over 0.93.

I. INTRODUCTION

Software-defined networking (SDN) is a paradigm that centralizes network control decisions and enables the network to be intelligently and centrally programmed. These characteristics simplify network management and provide tools for infrastructure sharing [START_REF] Mckeown | OpenFlow: Enabling Innovation in Campus Networks[END_REF].

SDN centralization provides advantages and disadvantages in terms of network security: on one side, the controller's global view has been used to develop new security strategies [START_REF] Ahmad | Security in Software Defined Networks: A Survey[END_REF], on the other side, the controller is a single point of failure, which turns SDN-based networks prone to DoS attacks [START_REF]New-flow based ddos attacks in sdn: Taxonomy, rationales, and research challenges[END_REF] [START_REF] Rawat | Software defined networking architecture, security and energy efficiency: A survey[END_REF]. In the case of resource constrained networks, as wireless sensor networks for Internet of things, SDN vulnerabilities are critical since there are less resources to detect and mitigate attacks. Consequently, current standard SDN security solutions adaptation is not trivial.

Since SDN centralizes the control logic of the network, most of works in the literature propose centralized security solutions. This has benefits, such as a global view of the network and high processing power, but it also requires a constant communication between the network devices and the controller. This means more energy and communication resources consumption. To address this issue, we propose a lightweight DoS attack detection algorithm using change point (CP) analysis to detect anomalies in the network behavior. We execute our proposal using a distributed approach, running the detection algorithm on individual resource constrained nodes, avoiding the packets overhead caused by the centralization.

We simulate grid topologies of 100 nodes, where 10% of nodes are attackers. Our main results show that individual nodes can detect a DoS attack and identify the attacker itself with a probability equal or over 0.93, when being close to the attacker. In addition, we investigate trade-offs between a fully decentralized and a hybrid approach.

II. RELATED WORK

In this section, we analyze recent works that propose DoS attacks detection solutions for SDN-based networks. We compare our proposal and the state mainly based on detection performance and resources limitations.

Machine learning is a popular approach used for security in SDN since the controller has access to traffic information that could be used to train the algorithms. Bhunia and Gurusamy [START_REF] Bhunia | Dynamic attack detection and mitigation in IoT using SDN[END_REF], Ravi and Shalinie [START_REF] Ravi | Learning-driven detection and mitigation of ddos attack in iot via sdn-cloud architecture[END_REF], and Jia et al. [START_REF] Jia | Flowguard: An intelligent edge defense mechanism against iot ddos attacks[END_REF] proposals have in common that all of them obtained high detection rate results, i.e., higher than 90%, using machine learning techniques. On the other hand, none of these three proposals considered resource constraints. The main reason is because these are OpenFlow-based or require high traffic of packets to monitor the network. Some proposals focused on resource limitations. Yin et al. [START_REF] Yin | A DDoS Attack Detection and Mitigation With Software-Defined Internet of Things Framework[END_REF], Miranda et al. [START_REF] Miranda | A collaborative security framework for software-defined wireless sensor networks[END_REF], and Wang et al. [START_REF] Wang | ETMRM: An Energyefficient Trust Management and Routing Mechanism for SD-WSNs[END_REF] proposed more lightweight security solutions, but at the cost of detection rates below 90%. However, these works proposed multiple types of attack detection and attacker identification algorithms.

Commonly, security in SDN is centralized, but nevertheless, there are distributed-based proposals in the literature. One reason to use distributed approaches is to avoid control overhead that could saturate vital control links. To address this shortcoming, Fawcett et al. [START_REF] Fawcett | Tennison: A distributed sdn framework for scalable network security[END_REF] proposed Tennison, a framework for scalable network security based on multi-level flow monitoring. Distributed approaches have been used also to detect anomalies in local sub-networks [START_REF] Shu | Collaborative intrusion detection for vanets: A deep learning-based distributed sdn approach[END_REF] The main shortcoming in the state of the art is the tradeoff between detection rate and solution complexity. The proposals that attained high detection rate were not suited for resource constrained networks, and proposals that considered resource limitations did not attain high detection rates. In this work we propose a DoS attack detection algorithm for SDN resource constrained networks, lightweight enough to run on individual resource constrained devices. Results show a detection probability comparable to centralized proposals, but reducing packets traffic, a key shortcoming in centralized solutions.

III. SDN SECURITY VULNERABILITIES

As explained in Section I, network control centralization and planes separation are fundamental enablers of SDN programmability. On the other hand, these traits turn the network vulnerable to denial of service (DoS) attacks,

In SDNs, the attackers can reach the control plane directly through the controller(s) or through network devices. Control packets flooding attacks are common since these packets have to reach the controller to be processed, which can lead to processing and communication resources exhaustion. The attackers are also able to mislead other network devices and induce them to flood the network.

The SDN controller needs topology information to operate. To this end, the network devices send neighborhood information to the controller for configuration and control decisions. In the case of wireless SDN networks, attackers may hear this information and use it to mislead the controller to take wrong routing decisions.

In the case of SDN resource constrained networks, these attacks target specific characteristics. Attackers can launch control plane attacks to saturate flow tables and buffers of devices with low storage capacity. A saturated node may not have space to forward new packets or receive new routing rules. This will trigger a series of packets retransmissions, which means more energy, processing and communication resources consumption. Since the network operation depends on the controller, if this do not take actions, the network devices can exhaust all their resources.

In a previous work [START_REF] Segura | Understanding the Performance of Software Defined Wireless Sensor Networks Under Denial of Service Attack[END_REF], we analyzed the impact of DoS attacks in SDWSNs. One of these attacks was the false data flow forwarding (FDFF) attack.The FDFF attack targets the controller via network's devices. First, the attacker sends data packets with unknown flow identifiers to its neighbors. The neighbors receive the packet and check the flow table to determine the action required, without success, thus they ask a rule to the controller by sending a flow rule request packet. The controller receives this packet, calculates the rule and replies sending a flow setup packet. This attack increases the control overhead and the processing overhead on the nodes in the neighborhood and on the controller. Also, after several repetitions, this attack can saturate the neighbors' flow tables.

IV. DISTRIBUTED DOS ATTACK DETECTION

From [START_REF] Gustavo | Denial of Service Attacks Detection in Software-Defined Wireless Sensor Networks[END_REF] and [START_REF] Segura | Multimetric online intrusion detection in software-defined wireless sensor networks[END_REF] we know that our CP detection algorithm, based on [START_REF] Skaperas | Early Video Content Popularity Detection with Change Point Analysis[END_REF] [START_REF] Skaperas | Real-Time Video Content Popularity Detection Based on Mean Change Point Analysis[END_REF], is able to detect FDFF attacks with a probability over 0.96, and identify the type of attack with a probability exceeding 0.89. In this work, we go further and execute the CP detector in a distributed approach, this means, running on individual nodes. Our objectives are: first, to evaluate the performance of the CP detector in resource constrained devices, and second, study the tradeoff when running the detectors on every node in the network and running it in clusters.

A. Change point detection

The problem formulation exploits recent results [START_REF] Skaperas | Early Video Content Popularity Detection with Change Point Analysis[END_REF], [START_REF] Skaperas | Real-Time Video Content Popularity Detection Based on Mean Change Point Analysis[END_REF] on non-parametric real-time CP detection. We adapted the hybrid offline-online proposal to an entirely online detector [START_REF] Segura | Centralized and distributed intrusion detection for resource constrained wireless sdn networks[END_REF].

To outline the online CP algorithm, let {X n : n ∈ N} be the time series of the metric monitored. Using Wold's theorem we can assume that, for X 1 , ..., X N , each sample is expressed as X n = µ n + Y n , where {µ n , n ∈ N} is the mean of the time series and {Y n : n ∈ N} is a random zero mean term, so that we can rewrite X n as:

X n = µ + Y n , n = 1, . . . , m + k * -1 µ + Y n + I, n = m + k * , . . . (1)
where µ, I ∈ R r , represent the mean parameters before and after the unknown time of possible change k * ∈ N * respectively. The term m denotes the length of an initial period assuming no change on the mean value, i.e,

µ 1 = • • • = µ m .
During this period, our detector "learns" in real-time the statistics of the observed time series, and, the mean value in particular. Finally, the statistical hypothesis test is articulated as: H 0 : I = 0, H 1 :

I = 0.
The online analysis is a stopping time stochastic process defined as:

τ (m) = min{l ∈ N : T S on (m, l) F (m, l)}, ∞, if T S on (m, l)< F (m, l) ∀l ∈ N, (2)
where T S on (m, l) is the detector, calculated online for every l, and F (m, l) is the given threshold; with properties lim m→∞ P r{τ (m) < ∞|H 0 } = α, ensuring that the probability of false alarm is asymptotically bounded by α ∈ (0, 1), and, lim m→∞ P r{τ (m) < ∞|H 1 } = 1, ensuring that under H 1 the asymptotic power is unity. Under these conditions, F (m, l) = cv on,α g(m, l), where the critical value cv on,α is determined from the asymptotic distribution of the detector under H 0 and the asymptotic behavior achieved by letting m → ∞. The weight function is defined as,

g(m, l) = √ m 1 + l m l l + m γ (3)
where the sensitivity parameter γ ∈ [0, 1/2).

The online algorithm uses the standard CUSUM detector [START_REF] Fremdt | Asymptotic distribution of the delay time in page's sequential procedure[END_REF], with test statistic denoted by T S ct on . Its corresponding critical value is denoted by cv ct on,α and the stopping rule by τ ct (m). The sequential CUSUM detector is denoted by E(m, l) = X m+1,m+l -X 1,m . The standard CUSUM test is expressed as:

T S ct on (m, l) = l Ω -1 2 m E(m, l), (4)
where Ω m is the estimated long-run covariance, defined as in [START_REF] Rawat | Software defined networking architecture, security and energy efficiency: A survey[END_REF], that captures the dependence between observations. Then, the stopping rule τ ct (m), is defined as:

τ ct (m) = min{l ∈ N : T S ct on (m, l) 1 ≥ cv ct on,α g(m, l)}, (5)
where the 1 norm is involved to modify T S ct on so that it can be compared to a one dimensional threshold function. The critical value, cv ct on,α , is derived from the asymptotic behavior of the stopping rule under H 0 :

lim m→∞ P r{τ (m) < ∞} (6) = lim m→∞ P r sup 1 l ∞ T S ct on (m, l) 1 g(m, l) > cv ct on,α = α (7)

B. DoS attack detection: implementation

In our distributed proposal, every node is able to monitor a time series using its own metrics and execute the CP detector algorithm, and also, each node is able to send metrics samples to a cluster head (CH). In the second case, the CH is in charge of constructing the time series of the cluster and execute the detection algorithm. Such approaches, in which hardware behavior is monitored to identify anomalies, hints to further integration with other approaches for the domain of physical layer security [START_REF] Chorti | Physical layer security: A paradigm shift in data confidentiality[END_REF] and introducing security controls at all layers.

A security application was programmed in every node. This application manages the sampling and the detection algorithm execution. The algorithm initiates constructing a time series of 200 samples (m = 200), from where it extracts the statistical information that will use during the online CP detection. When the online phase starts, the algorithm continues storing samples. In the case no CP is detected during the first 50 samples (l = 50), these samples are added to the 200 samples taken before to extract new statistical information. This, process is iterated every 50 samples, but, if a CP is detected, the application rises an alarm and informs the controller about the situation.

Since the FDFF attack increases the message exchange activity, we decided to monitor the transmitting time on every node, i.e, the number of ticks the radio module remained turned on transmitting packets. From previous results [START_REF] Segura | Multimetric online intrusion detection in software-defined wireless sensor networks[END_REF], we learned that γ = 0 and m = 200 maximize the detection rate. In this work we used these values as well.

The security module was implemented in C language using Contiki-3.0 [START_REF] Dunkels | Contiki-a lightweight and flexible operating system for tiny networked sensors[END_REF], an operating system for WSN and IoT, and IT-SDN [START_REF] Alves | The Cost of Software-Defining Things: A Scalability Study of Software-Defined Sensor Networks[END_REF], an SDWSN framework developed by our reserch group. For transmitting time sampling

V. RESULTS AND ANALYSIS

Our analysis follows two approaches: detection performance and implementation overhead. For detection performance we analyzed the attack detection probability and attacker identification probability. For implementation, we analyzed the packets overhead and memory usage. Both scenarios were simulated on Cooja [START_REF] Osterlind | Cross-Level Sensor Network Simulation with COOJA[END_REF], simulating grid topologies of 100 nodes, where 10% of nodes are attackers and emulating sky motes.

A. Detection performance

Fig. 1 shows the detection probability heatmap for both the distributed and hybrid approaches. For the case where each individual node is running the detector (Fig. 1a) the attackers' position is represented with an "A". For the case where the detector is running on clusters (Fig. 1b), the heatmap shows also the number of attackers in each cluster.

In Fig. 1a we observed that the detection probability is higher around the attackers and around the controller, with values between 0.93 and 1.00. This is the behavior we expected since the FDFF attack targets the control plane through the attackers' neighbors, thus the attack has a direct impact on the transmitting time mean value of these nodes. The CPs detected on the nodes at two or more hops from the attackers reflect the impact of this attack in the whole network. The increase in the transmitting time of these nodes is caused mainly by the increase in the control packets forwarding.

The detection probability results for the clustering case in Fig. 1b showed that in the clusters where there are one or two attackers, the detection probability was equal or above 0.97. In the clusters without attackers, the detection probability was between 0.63 and 1.00. Similar to the case running the detector on individual nodes, high detection probability results in clusters without attackers mean that this attack impacts on the whole network.

The detection probability results provided two important insights: (i) we were able to detect an FDFF attack when monitoring the transmitting time in either individual nodes or clusters, and (ii) the detection probability was lower on the nodes at two or more hops from the attackers. Based on (ii), we implemented an attacker identification algorithm for the case running on individual nodes. First, every time a node receives a data packet with an unknown flow, it saves the identification address of the sender in a vector for suspects. Then, if a node detects a CP on the transmitting time, it sends an alarm to the centralized security module and informs the address of a suspect. To determine the suspect, the node checks the last ten addresses saved in the vector and choose the one with the highest frequency. We chose ten samples because, according to [START_REF] Segura | Multimetric online intrusion detection in software-defined wireless sensor networks[END_REF], the slower detection when γ = 0 takes around ten samples. When the security module receives the alarm, sets a counter for every suspect. If the counter is equal to the number of neighbors of the corresponding suspect, the suspect is declared as attacker.

The heatmap in Fig. 2 showed that using this algorithm we were able to correctly identified all the attackers with a probability equal or above 0.93. In addition, the false positives were equal to zero. On the other hand, this was possible only running the detector on every individual node. In the clustering approach, groups without attackers inside also obtained high detection rates, which excluded the possibility of tracking the attacker based only on the alarm received from the cluster.

B. Implementation comparison

In the clustering approach, the cluster head is in charge of constructing the time series for the whole cluster and execute the CP detector. To accomplish this, all nodes have to monitor their transmitting time and send a sample to the cluster head, periodically. The cluster head sums up all the samples received per period and the result is a new sample of the cluster's time series. Conversely, when the CP detector is running on individual nodes, each one constructs its own transmitting time time series and executes the CP detector.

The clustering approach has the benefit that it reduces the memory and the processing overhead of the network, since the detector is running on the cluster heads only. On the other hand, it increases the packets traffic, since all nodes now have to periodically send a transmitting time Fig. 2: Attacker identification probability: sample to the cluster head. In terms of memory, the time series construction and the CP detector implementation, for one metric, represents 5956 B. In our specific case using sky motes, this value represents 12.40% of the total memory.

One security risk when using clustering approaches, is the possibility of the cluster head to be an attacker. One way to solve this is using secure cluster head selection algorithms [START_REF] Shankar | Hybrid model for security-aware cluster head selection in wireless sensor networks[END_REF], but this requires more memory, processing, and communication resources, which are already scarce in our case. One benefit of running the detector on every node is that we avoid this risk since the attack detection does not depend on one or few nodes only.

VI. CONCLUSION

In this work we propose a distributed DoS attack detection for software-defined resource constrained wireless networks, based on CP detection. Our proposal is lightweight enough to run on very limited devices and detect FDFF attacks with a probability between 0.93 and 1.00, comparable to detection results in centralized proposals.

This proposal was evaluated running the detector on every node and running the detector in clusters. Results showed that both approaches obtained high detection rate. Additionally, we were able to identify the attackers without increasing the packets traffic when running the detector on every node. Comparing the implementations of both approaches, the clustering approach requires less memory on every node, while running the detector on every node reduced the packets traffic.

As future work, we envisage to develop a full implementation of distributed and centralized approaches and compare their performance based on security and network performance metrics.

Fig. 1 :

 1 Fig. 1: Distributed FDFF attack detection

ACKNOWLEDGMENTS

This study was financed in part by the Coordenac ¸ão de Aperfeic ¸oamento de Pessoal de Nível Superior -Brasil (CAPES) -Finance Code 001 and by the ELIOT project (ANR-18-CE40-0030 / FAPESP 2018/12579-7). Gustavo A. Nunez Segura is supported by Universidad de Costa Rica.