Matteo Sensi

Carole Baffert

Vincent Fourmond

Luca De Gioia

Luca Bertini

Christophe Léger

Photochemistry and photoinhibition of the H-cluster of FeFe hydrogenases

come

Introduction

Dihydrogen is a very common energy vector in Nature, as most microorganisms either oxidize or produce H 2 in relation to their energy metabolism. These reactions are catalyzed by enzymes called hydrogenases 1 and nitrogenases.

H 2 oxidation occurs in relation to bacterial respiration: it is coupled to the reduction of a number of electron acceptors such as sulfate in sulfate reducing bacteria or O 2 in Knallgas bacteria, and the free energy of the reaction is used for producing chemical energy in the form of ATP.

The dihydrogen that is consumed by this process can be produced in a number of ways, either in the same microorganism or in a microorganism from the same consortium.

The most prominent H 2 -producing metabolic pathway is dark-fermentation 2 , as e.g. occurs in strictly anaerobic Clostridial species [START_REF] Latifi | [END_REF] or in enteric bacteria such as Escherichia coli 2 .

H 2 is also a byproduct of the reduction of N 2 to NH 3 , a key step in the global nitrogen cycle which occurs e.g. in cyanobacteria: the latter reaction is catalyzed by the enzyme nitrogenase, which evolves at least one molecule of H 2 per N 2 reduced 4,5 .

Biological splitting of water into O 2 and H 2 can also be observed in photosynthetic microorganisms such as green algae and cyanobacteria 6 , but the combination of O 2 evolution and H 2 production is not a natural physiological pathway 7 , and it is severely limited by the oxygen sensitivity of the enzymes involved in H 2 production: hydrogenase or nitrogenase 6,[8][9][10][11] . Efforts have been made to genetically engineer the enzymes [12][13][14] or the strains 10 of these photosynthetic microorganisms, aiming for continuous, light-driven biological H 2 production [9][10][11][12][13][14][15] . Another strategy to work around the O 2 sensitivity problem in microalgae consists in alternating aerobic and H 2 production phases by regulating the activity of Photosystem II; reducing equivalents are stored in the form of starch during the aerobic phase and then fed into hydrogenase via Photosystem I 11,16 . Alternatively, spatial separation of water oxidation and H 2 production can also prevent nitrogenase inhibition by O 2 in cyanobacteria 11 .

Biological, photosynthetic H 2 -production is therefore not easily achieved 17 , but chemists have long worked towards the goal of water-splitting-based hydrogen production by designing original semi-artificial systems that interface hydrogenases or hydrogenase mimics with either photoactive materials or molecular photosensitizers, including the photosystems involved in natural photosynthesis [18][19][20][21][22][23][24][25][26] .

Regarding the hybrid devices that include a photosynthetic enzyme, limitations may arise from the photodamage of the photosystems, which has been thoroughly studied by the photosynthesis community [27][28][29][30] . In vivo, the effect of the photodamage is balanced by de novo synthesis and repair mechanisms, but in vitro the damage decreases the stability of the hybrid systems [START_REF] Plumeré | Biophotoelectrochemistry: From Bioelectrochemistry to Biophotovoltaics[END_REF][START_REF] Zhao | [END_REF] .

An Achilles' heel of the artificial photoproduction systems that has received much less attention is the electron acceptor side: there is published evidence that some hydrogenases are sensitive to light irradiation, but this was rarely considered as a potential weakness 25,33,34 . This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D Nature makes three distinct enzymes to heterolytically split molecular H 2 . The so-called "NiFe hydrogenases" and "FeFe hydrogenases" both catalyze H 2 oxidation and production at a binuclear metal center and have been considered in relation to the design of artificial photosynthesis systems. In contrast, "Fe hydrogenases" catalyze hydride transfer from H 2 to methenyl-tetrahydromethanopterin, an organic substrate involved in microbial methane production. All three types of active sites contain at least one Fe-carbonyl bond 33 , and are sensitive to light but in a different manner, and according to mechanisms that have not been fully identified.

The active site of Fe hydrogenases 34 is inactivated by exposure to UV-A (320-400 nm) or blue-light (400-500 nm), as a result of Fe release from the cofactor 35 . The enzyme is protected from light-induced inactivation by the binding of the inhibitor CO 35 and, in the particular case of the enzyme from Methanothermobacter marburgensis, by its oligomerization 35 . A key aspartate residue and some protein rearrangements are involved in the Fe hydrogenase protection from light-induced decomposition of the active site 36 .

The active site of NiFe hydrogenases exists in the dark in a number of forms, some of which (the so-called Ni u -A, Ni r -B, Ni-C and Ni a -L states) are EPR active. Ni u -A and Ni r -B correspond to oxidized, inactive states of the enzyme that disappear upon reductive activation of the enzyme. The Ni u -A 37 and Ni a -C [38][39][40][41] states and the EPR-silent Ni a -S 42,43 and Ni-S-CO [44][45][46][47] states are all altered upon irradiation with visible light under cryogenic conditions. Ni a -C is an intermediate in the catalytic cycle, which bears a bridging hydride. Ni a -C is photoconverted under cryogenic temperatures into Ni a -L, where the electrons are formally stored on the nickel ion. Additionally, this photoconversion results in the change of the protonation state of a terminal Ni-coordinated Cys residue 48 . The photoproduct Ni a -L can also be formed in the dark at ambient temperature 49 . The CO-inhibited Ni a -S state is light-sensitive: the extrinsic Ni-bound CO dissociates under illumination at cryogenic temperature (<120 K) restoring the catalytically relevant Ni a -S intermediate [44][45][46][47] . Illumination with blue light accelerates the reactivation of the Ni r -B state of the enzyme from Aquifex aeolicus at room temperature 50 . A mechanism for the removal of the bridging OH ligand from the N i r-S inactive ready state upon irradiation has been suggested 51 . Re-activation can also be triggered by light under cryogenic conditions in the protein crystals 52 . We are not aware of any evidence that NiFe hydrogenases are damaged upon irradiation.

In the context of designing H 2 photoproduction devices, FeFe hydrogenase and their synthetic mimics are preferred because they are deemed more active for H 2 production 53,54 .

In this review we focus on FeFe hydrogenases and their biomimetic models, and we discuss the various effects that have been observed upon irradiation.

The structure of FeFe hydrogenases

The photochemical properties of FeFe hydrogenases are determined by the electronic structure of their inorganic cofactors. The active site of FeFe hydrogenases 55 , called the H-cluster, consists of two subclusters bridged and attached to the protein matrix by a cysteine residue: a canonical [4Fe4S] H cluster and a [2Fe] H site; the two Fe ions are called proximal, Fe P , and distal, Fe D , in reference to the distance to the [4Fe4S] H subcluster. Each Investigations on the photochemical properties of the H-cluster aimed at understanding the effect of light irradiation and the nature and reactivity of the photoproducts. Fast or ultrafast processes are difficult to observe, but the dynamics of the excited states can be studied using quantum chemistry with time-dependent DFT (TDDFT).

In general, electronic excited states are always calculated from the ground state wave-function, whose characterization is required before the effect of excitation can be investigated. The ground state electronic structure of the H-cluster in the H ox and H ox -CO states can be described as a three-spin system: the [2Fe] H site has a spin of 1/2, each [2Fe2S] layer of the [4Fe4S] H cluster has S = S(2)+S(5/2) = 9/2 (figure 2). The iron atoms of the binuclear subcluster, in both H ox and H ox -CO states of the H-cluster, are in the Fe I Fe II redox configuration. The attribution of the redox state to Fe P and Fe D is not trivial: FTIR 66,67 , This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D EPR 77 and DFT 78 studies suggest that Fe D is the Fe I while ENDOR, HYSCORE 79 and FTIR spectroelectrochemistry 61,65,80 data support the paramagnetic Fe I nature of Fe P . The [4Fe4S] H 2+ cluster consists of two [2Fe2S] layers coupled antiferromagnetically, each formed by one Fe II and one Fe III high-spin ions, coupled ferromagnetically, giving an overall low-spin state 58,78,81,82 .

Figure 2 -Three-spin model of the H-cluster in the H ox state.

The frontier molecular orbitals (FMOs) of the H-cluster in the ground state (GS) obtained from DFT calculation are localized on the [4Fe4S] H cluster and/or the [2Fe] H site, depending on the state of the H-cluster 81,83 We have calculated the frontier molecular orbitals (FMOs) of H ox -CO in the ground state 84 . In this case, single occupied HOMO and LUMO are localized on the [4Fe4S] H and on the [2Fe] H respectively. We described the H-cluster as a dyad system that consists of two fragments covalently bound to one another: the [4Fe4S] H cluster and the [2Fe] H site. As shown in scheme 1, the excited states decay of this system can be described in terms of a charge-transfer (CT) or a localized mono-electronic transition depending on the orbital composition of the FMOs involved 85,[START_REF] Chongdar | This is a post-peer-review, pre-copyedit version of an article[END_REF] .

Scheme 1 -The electronic structure of the FeFe hydrogenase H-cluster, seen as a molecular dyad. The scheme of the excitations is proposed on the basis of the frontier MO composition of the H ox -CO model described in ref. 84 . The lower part is adapted with permission from ref 84 , copyright 2016 American Chemical Society.

The above-described H-cluster defines the large family of FeFe hydrogenases. These enzymes are present in many microorganisms. They differ in terms of protein sequences, quaternary structures and cofactor content [START_REF] Poudel | [END_REF]88 . The immediate environment of the H-cluster is mostly conserved, although some variations have been observed and may induce unusual catalytic properties 55,85,[START_REF] Chongdar | This is a post-peer-review, pre-copyedit version of an article[END_REF] . Non conserved residues that are remote from the active site may This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D also affect active site chemistry, as demonstrated using site-directed mutagenesis in the particular case of Clostridium Beijerincki FeFe hydrogenase 89 . They embed a variable number of accessory FeS clusters which are used to mediate long range intramolecular electron transfer: there are 2 accessory clusters in the FeFe hydrogenases from D. desulfuricans (DdHydAB), D. vulgaris (DvHydAB) and M. elsdenii (MeHydA) FeFe hydrogenases, 4 in the FeFe hydrogenases from Clostridium acetobutylicum (CaI), and Clostridium pasteurianum (CpI), and none in the CrHydA1 enzyme from Chlamydomonas reinhardtii (figure 1).

The effect of light on the CO-inhibited state of FeFe hydrogenases

The reactivation of the CO-inhibited FeFe hydrogenase by light has been studied since the 1930s, initially to learn about the metal content and structure of the active site and the binding site of hydrogen. More recent investigations of the effect of irradiation have focused on understanding the properties and the reactivity of the active site (Table 1).

Early functional studies of the effect of light on CO-inhibited hydrogenases

In 1933, Kempner and Kubowitz showed that the production of H 2 by Clostridium butyricum is inhibited by CO and illumination reverses the inhibition 90 . At that time neither the structure of the enzyme nor the presence of nickel and/or iron in the active site had been established; the flash photolysis of CO had been used to study the reaction of myoglobin with small ligands 91,92 and the light-dependent CO-inhibition of hydrogenases was taken as an indication of the presence of an iron atom at the active site, as in myoglobin.

The reversal of CO inhibition by white light was also observed with intact cells and cell-free extracts of other bacteria [93][94][95][96][97] . Some of these early observations are ambiguous because some of the bacteria that were investigated, e.g. D. desulfuricans, produce FeFe and NiFe hydrogenases, both of which are light-sensitive in the CO-inhibited state. The effect of irradiation was probably also modulated by the permeability of the biological membranes.

Rudolf K. Thauer and coworkers performed the first kinetic study of FeFe hydrogenases inhibition by CO in vitro and also observed the reversal of the inhibition by light 98 . The authors qualitatively observed that after exposure to CO, the H 2 oxidation activity of the enzyme is ten times greater under the light of a halogen lamp than in the dark, which was interpreted as evidence of the photolysis of the CO-enzyme complex. The effect of light was fully reversible. The rates of binding and release of CO under illumination were not measured.

Spectroscopic and crystallographic investigations

EPR and FTIR spectroscopic techniques have later been used to show that illumination affects the structure of the H-cluster in the H ox -CO state. Under cryogenic conditions, one CO ligand is photo-dissociated, but the structure of the photoproduct depends on temperature. At room temperature the illumination allows the exchange between intrinsic and extrinsic CO, demonstrating the dynamics of the [2Fe] H site.

This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy

& Fuels (2021) 10.1039/D1SE00756D

Low temperature

Peters and Lemon obtained the first crystallographic structure of the FeFe hydrogenase of Clostridium pasteurianum (CpI) 99 , and in 2000 they reported the effect of irradiating the crystal of the CO-inhibited state with a red laser (635.8 nm) under a liquid nitrogen stream 100 . The exogenous CO binds on the axial site of Fe D and partial dissociation is observed after irradiation. The difference between the electron density maps evidences an increased distance between the S of a conserved methionine residue and the bridging CO, resulting from an increased electron density on the bridging CO after the photolysis of the apical CO 100 .

In the late 1980's, Patil 101 , Kowal 102 and coworkers had used EPR to examine the effect of light on the enzymes from D. vulgaris (DvHydAB) and C. pasteurianum (CpI), respectively, under cryogenic conditions. Patil and colleagues observed that the H ox -CO EPR signal at g = 2.06 signal is converted to the rhombic g = 2.10 signal, assigned to H ox , upon irradiation with a 200 W mercury arc lamp at 9 K. Heating the sample to 200 K gives back the H ox -CO signal, showing that the CO ligand remains in the environment of the H-cluster and rebinds if the temperature is high enough that the corresponding energy barrier can be overcome. According to Kowal et al., irradiation at 8 K produces H ox , but irradiation at 30 K produces a signal at g = 2.26, 2.12 and 1.89, falsely suggesting the presence of a low spin Fe III center, at a time when the structure of the H-cluster was unknown.

Changes in H-cluster structure have also also been detected using FTIR spectroscopy. Peters, Bagley and coworkers observed that the irradiation of CpI H ox -CO with a 300 W xenon lamp at 6-14 K produces the H ox state (which decays back to H ox -CO above 150 K), but a second species is produced upon irradiation above 20 K; it is characterized by the loss of the bridging CO band at 1809 cm -1 and decays above 80 K by CO rebinding 67 . The authors hypothesize that the later species gives the g = 2.26 rhombic signal in the experiments of Kowal et al. 101,102 , and that the photodissociation involves the bridging or apical CO depending on temperature (see fig. 3). Working with D. desulfuricans hydrogenase (DdHydAB), Roseboom et al. also observed by FTIR that illuminating H ox -CO at 15-30 K with an halogen lamp or with a more intense light (xenon lamp), causes the disappearance of both the bands of extrinsic CO and of the bridging CO, or only of the bridging CO in molecules where the extrinsic CO is already photo-dissociated 66 . The removal of the bridging CO from H ox -CO at cryogenic temperatures has no effect on the frequencies of the other intrinsic COs and on the charge densities of the irons, again in contrast with the results in ref 67 . It is unclear if the observed differences result from the enzymes from C. pasteurianum and D. desulfuricans behaving differently, or the different irradiation times (10 minutes versus 20 minutes to hours) or light power (300 W vs 150 W) or lamp emission spectra. Based on an EPR and ENDOR investigation, Lubitz and coworkers concluded that the "second" photoproduct (H LI ) is characterized by an unpaired spin density localized on the diiron site, in particular over the CN -ligands, resulting from a structural rearrangement which cannot be merely the dissociation of bridging CO 103 . In all cases, all the signals obtained after irradiation decay upon warming the sample in the dark at 150 K. 67 and Roseboom 66 . Irradiation of H ox -CO state at temperatures between 30 and 150 K produces essentially H ox , while no CO photolysis is observed at T>150 K due to rapid recombination of the CO ligand. Additionally, all light-induced changes can be completely reversed by warming the sample to 150 K in the dark.

Room temperature

In experiments carried at room temperature, Albracht et al. observed that the illumination with an halogen lamp has no effect on the EPR and FTIR spectra of the H ox -CO and H inact states, but in samples prepared with 13 CO illumination shuffles the two distal CO ligands bound to Fe D 59,66 .

Stripp et al. used this effect to study the dynamics of the [2Fe] H site by isotope exchange with real-time attenuated total reflection FTIR (ATR-FTIR) 106 . The exposure of CrHydA1 in the H ox -CO state to exogenous 13 CO results in a shift of the IR bands that is ascribed to the exchange of one of the intrinsic CO ligands on Fe D with the labelled ligand. Subsequent irradiation with red light (640 nm) and then blue light (460 nm) induces the formation of a cluster binding four labeled 13 ligand in the apical position, instead of the consensual equatorial position. In contrast, other spectroscopic 107 and DFT studies [108][109][110] have considered unlikely the coordination of the distal CN -on the apical binding site.

Kinetics

The kinetics of CO binding or release has been examined using time-resolved IR spectroscopy, to focus on rebinding of CO after photodissociation, and also under turnover conditions by protein film voltammetry, to examine how irradiation affects the kinetics of inhibition by CO.

Regarding the rebinding experiments, Mirmohades et al. 111 irradiated a C. reinhardtii

FeFe hydrogenase 1 sample in the H ox -CO state by flashing a UV laser (10 times 10 ns flash, λ = 355 nm), at room temperature, and used time-resolved nanosecond mid-infrared spectroscopy to detect CO photo-dissociation and rebinding. The quantum yield of conversion is around 10% for this enzyme, and around 8% of the total sample was converted into H ox after flashing. The kinetics of rebinding was not 1st-order, suggesting that the dissociated CO may reach several binding sites in the protein before rebinding to Fe D , and the half-life of the rebinding of CO was 13 ± 5 ms.

In protein film voltammetry, the enzyme is adsorbed onto an electrode, and the electrons produced upon H 2 oxidation are directly transferred to the electrode. The resulting current is proportional to the turnover frequency, and the variations in current report on any changes in activity [112][113][114] which may result from CO binding or release. A light source can be directed towards the electrode that supports the enzyme, to examine how irradiation affects the activity or the inhibition kinetics 112,113 .

Parkin 115 , Goldet 116 and coworkers in the Armstrong group have first used this approach to qualitatively evidence the effect of light on the recovery of activity of the CO-inhibited enzyme. The illumination of the enzyme with the white light from a halogen lamp (10 W and 20 mW/cm 2 in ref 115 , 150 W in ref 116 ) enhances the recovery of activity, but only under oxidizing conditions; no effect is seen at low electrode potential, when the enzyme is evolving H 2 . No difference was noted between C. acetobutylicum HydA1 and C. reinhardtii HydA1. The rate of the process was not determined.

We also used electrochemistry to examine how light affects the CO binding and release kinetics, aiming for a quantitative description of the photoreactivity. We irradiated the FeFe hydrogenases from M. elsdenii, C. acetobutylicum and C. reinhardtii (HydA1) adsorbed onto electrodes with monochromatic lasers, varying the wavelength (405, 450, 532 and 635 nm) and the incident power density (from 0 to 60 μW/cm 2 ). We obtained the values of the CO binding ("in'') and release ("out") rate constants using previously developed methods 14,[117][118][119] . We found that irradiation has no effect on the CO binding rate, but blue or violet light efficiently dissociates the exogenous CO. From the slope of the linear change in k out against light power, we estimated the value of the molecular absorption cross-section relative to the photodissociation of CO, which appears to be similar to the values deduced from the extinction coefficient of H-cluster biomimetic complexes; this shows that most of the absorbed blue-violet light leads to the dissociation of CO.

We combined the experimental results with the conclusions from TDDFT calculations to obtain an electronic and mechanistic description of the process. The CO inhibited H ox -CO This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D form of the H-cluster is an oversaturated organometallic complex since one of the two Fe atoms formally counts 19 valence electrons. Our TDDFT calculations 120 show that H ox -CO absorbs at a lower wavelength than H ox , and evidence a resulting CO photolysis pathway. Indeed, starting from the electronic structure of the H-cluster 83 , the H ox -CO calculated spectra showed that only the Fe 4 S 4 →Fe 2 S 2 charge-transfer and Fe 2 S 2 localized excitations in the visible light range are active toward CO dissociation, in agreement with the experimental results. The exploration of the potential energy surfaces of the excited states as a function of the Fe-CO bond length clearly evidenced a non-radiative pathway characterized by a series of internal conversions (IC) that end with the release of the CO in the apical position. Using TDDFT geometry optimizations we also simulated one of these possible pathways for the dissociation of the apical CO 84 .

Photoinhibition

Early FTIR and EPR measurements showed DdHydABis damaged by the sunlight or by the white light from a Xe lamp, but not a halogen lamp. In contrast, recent PFV and FTIR studies of D. desulfuricans hydrogenase showed that inactivation occurs upon irradiation with a halogen lamp under oxidizing turnover conditions. Regarding CaI and CrHydA1 hydrogenases, no inactivation was observed upon illumination with a halogen lamp or lasers in the visible range, but UV-B irreversibly inactivates the enzyme under both oxidizing and reducing conditions (Table 1).

Figure 4 -Emission spectra of sunlight (red), a xenon lamp (black) and a halogen lamp (blue). Reprinted with permission from ref 121 , Copyright 2017 American Chemical Society.

Figure 4 shows the xenon and halogen lamps emission spectra. The halogen lamp spectra has an overall higher emission in the visible range above 600 nm than the xenon lamp but the latter has some emission peaks in the near-infrared, and higher UV components than the halogen lamp. The power density actually received by the sample also certainly matters 84 , since saturation of light absorption only occurs when the absorption cross section is very large; however, this information is rarely available in research papers.

The first evidence of the detrimental effect of light on FeFe hydrogenase was published by Albracht and coworkers in 2006 59 , in an IR and EPR investigation of D. desulfuricans hydrogenase (HydAB). At room temperature and upon irradiation with a UV-filtered xenon lamp, the EPR and FTIR 66 contributions from the H ox -CO state, which is a signature of H-cluster degradation: CO ligands are released from some degraded H-cluster and bind to other, still-intact enzyme molecules. This "cannibalization" process is also observed when the enzyme is exposed to O 2 122 . DdHydAB in the H red state was less sensitive to light. Under cryogenic conditions, illumination of the D. desulfuricans hydrogenase samples in the H ox -CO state with a xenon lamp leads to the disappearance of both the extrinsic and bridging CO bands from the FTIR spectra, whereas only the release of the extrinsic CO is observed using a halogen lamp. The specific effect of the two lamps suggests that the differences between the emission spectra of the halogen and xenon lamp are crucial. Xenon lamps have a UV component, which is almost completely absent in halogen lamps, but the authors used a "glass light guide" that eliminated "most of the UV contribution" 66 . We can only speculate that the different effects of the lamps are due to the higher intensity of the emission of the xenon lamp in the visible range or to a small fraction of not filtered UV light.

In 2016, Rüdiger and coworkers 123 also observed the light sensitivity of D. desulfuricans hydrogenase at room temperature, using protein film voltammetry and FTIR. The FTIR spectra of D. desulfuricans HydAB after irradiation presented the typical H ox -CO bands previously observed by Albracht et al. 59,66 In electrochemistry experiments, the authors illuminated the protein film with a 50 W halogen lamp positioned under the electrochemical glass cell, which we expect should block UV light. Under oxidizing conditions, constant illumination induces an inactivation followed by a partial recovery of the activity; another partial reactivation is observed when the light source is switched off. No such effect is observed under reducing conditions (figure 2b in ref 123 ). The authors propose that the photoinhibition of the enzyme on the electrode is due photodamage and cannibalization, and the subsequent reactivation would result from the dissociation of CO. However, and as pointed out by the authors themselves, the experimental conditions do not allow the small amount of CO resulting from enzyme degradation to bind to nearby enzyme molecules on the electrode surface. The two-step reactivation, under illumination and after light removal, is not explained.

In the above described experiments performed with D. desulfuricans hydrogenase, the role of the emission spectra of the lamp is not clear. Rüdiger et al. observed the inactivation of DdHydAB using a halogen lamp, whereas Albracht et al. observed the photodamage of the same enzyme upon irradiation with a xenon lamp but not with a halogen lamp. In independent electrochemistry experiments, Parkin et al. also saw no photoinhibition of D. desulfuricans hydrogenase by irradiating with a halogen lamp 115

(whose power was lower than that used by Rüdiger et al.).

Regarding the photoinhibition of the enzymes from C. acetobutylicum (HydA1) and C. reinhardtii (HydA1), Goldet and collaborators detected no photoinhibition upon irradiation with a halogen lamp; they only observed that irradiation slightly increases magnitude of the current and do not comment on whether or not it is possible that an increase in activity would result from irradiation heating the graphite electrode 116 . We consider it unlikely that irradiation would significantly perturb the catalytic cycle of hydrogenase: because of the small value of the absorption cross section of the enzyme in the visible range, photon absorption is a very rare event, unless a high power laser is used. On our side, we also used direct electrochemistry to look for evidence of photoinhibition of CaI and CrHydA1. We saw no detrimental effect of laser irradiation at 405, 450, 532 or 635 nm, but using a xenon lamp and various filters, we concluded that UV-B slowly and irreversibly inactivates the enzymes; this occurs under both H 2 oxidation and H + reduction conditions when the enzyme is adsorbed onto an electrode, and also when the enzyme is stored in solution. involving the intrinsic CO ligands that can be populated upon absorption in the UV domain.

This is a post-peer-review, pre-copyedit version of an article published in

Exploring the H ox excited states potential energy surfaces (PES) as a function of the Fe-CO and Fe-S bonds lengths of the [2Fe] H subcluster, we observed in silico the dissociation of the terminal CO ligands. The most stable photo-product has an all-terminal geometry without the bridging CO and with two equatorial COs, irrespective of the redox state of the initial species. The reduced form of this photo-product could easily form a stable dead-end with a bridging hydride, which we considered the likely structure of the inactive photoproduct 121 .

The reason why DdHydAB is apparently sensitive to white light, unlike homologous hydrogenases, is unknown 121 .

Photochemistry of mimics of the H-cluster

The hydrogenases scientific community has invested enormous efforts in the search of structural and functional analogues of the H-cluster. Many models have been synthetized, considering the diiron [Fe 2 (CO) 6 (pdt)] model as the parent 124,125,126 . Photochemical activation is a common tool for ligand substitution in this family of complexes 124 . Figure 5 shows three general mechanisms by which Fe 2 S 2 models can produce H 2 . A photocatalytic diad system is obtained when this cluster is conjugated to a photosensitizer (PS), such as a classical Ru(II) or Pt(II) complex with bis-or tris-dentate pyridine or porphyrins ligands (Fig 5b). 25,127 H 2 can then be produced in the presence of a less reducing electron donor, since the free energy for H 2 evolution comes in part from an electron or energy transfer process from PS to the cluster 128 .

An alternative approach consists in using a protonated Fe 2 S 2 model that can itself absorb light energy (Fig 5c). The advantage is to avoid the sensitization and all issues related to the chemical PS-model conjugation and the quantum yield of the PS-model quenching. The first photocatalyst of this type was the protonated Fe 2 (pdt)(μ-H)(CO) 4 (dppv) (where dppv= cis-1,2-C 2 H 2 (PPh 2 ) 2 ) investigated by Rauchfuss et al. [129][130][START_REF] Frederix | This is a post-peer-review, pre-copyedit version of an article[END_REF] In this case, light absorption induces a loss of the bridging character of the hydride, which is reminiscent of the behavior of the μ-CO in the CO-inhibited form of the H-cluster [129][130][START_REF] Frederix | This is a post-peer-review, pre-copyedit version of an article[END_REF] . The formation of a terminal hydride is favored by the presence of electron-rich phosphine ligands 132 . (Artero et al. 2015) .

The main focus of the photochemical studies was on the ultrafast (< 1 ns) photodynamics of the low-energy excitation, which is mainly dominated by CO photolysis, as is usual in metal-carbonyl complexes. The aim was to understand the effect of the ligands on the electronic properties of the model in relation to the nature of the photoproducts. The photostability of the models is crucial, if the goal is to wire them to photosensitizers for H 2 photoproduction. Flash photolysis kinetics, ultrafast pump/probe time-resolved IR (TRIR) spectroscopy and quantum chemistry calculations have been used to characterize the non-radiative transitions that occur upon UV-Vis light absorption.

The TRIR investigation of Fe 2 S 2 18-electron models evidenced at least one ultrafast CO dissociation process, with the formation of an 16-electron intermediate species, whose life time depends on the chemical nature of the ligands (e.g. 25 ps in the case of a phen substituted model 152 and much longer for a cyano-functionalized model 151 ). TDDFT calculations 144 suggests that dissociation after excitation in the UV domain results from the population of a dissociative Fe→CO charge transfer (CT) state, which involves an equatorial CO ligand.

When Fe 2 (CO) 4 (L) 2 (bridging ligand) models are considered, ligand isomerization is observed, which is due to the concomitant CO dissociation and Fe-S bond breaking during the non-radiative decay of the first excited state to the ground state 149,150 . The latter result is reminiscent of the photo-inhibition of the enzyme, where photochemical channels are populated when the enzyme is irradiated in the UV range. After CO photolysis, the formation of a long lived solvent adduct is observed 141 , as in the case of the mono cyano substituted models 154 , and geminate recombination gives back the initial model compound.

The excited state properties of the CO inhibited models have been investigated only theoretically, because over-saturated 19 electron complexes are not stable enough to be synthesized (see however ref. 158 ). According to TDDFT calculations, the Fe 2 (CO) 6 (μ-CO)(pdt) complex, a model of the H ox -CO form of the H-cluster, photo-releases a terminal CO after excitation in the visible domain in a lower energy range than observed with H ox models 153 .

Conclusions

The investigation of the photochemistry of the H-cluster biomimetic models has been important not only to design a better and efficient system for hydrogen photo-production, but also to shed light on the functionality of the active site itself. The spectroscopic, electrochemical and DFT studies on the photochemistry of FeFe hydrogenases contributed This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D not only to improve the development of hydrogen photo-production but also to obtain a fine description of the H-cluster electronic structure.

Our survey of H-cluster biomimetic model photochemistry has revealed that most of the complexes are photosensitive in the UV domain, as are the enzymes, and undergo ultrafast photodissociation of a terminal CO ligand. Following this, the unsaturated models form a solvent adduct that slowly decays to the initial model. The scrambling of the ligands is often observed along the vibrational cooling pathway, this process being favored by the transient vacancy of the CO ligand and the weakening of a Fe-S bond. The main difference between the photodynamics of the H-cluster and of its models comes from the ability of the latter to form a solvent adduct, whereas the H-cluster is shielded from the solvent.

The investigation on CO inhibited models has shown that the presence of an exogenous CO ligand lowers the energy range of the MLCT bands and of the CO photolysis channels. For this reason it appears that the ligand CO protects the [2Fe] H cluster since it is able to capture visible light, after which the system decays through photo-dissociative channels that involve a CO ligand, thus preventing irreversible photo-damage.

The literature on the photochemistry of the enzymes is related to photo processes that occur on longer time scales, and is relevant to the practical use of hydrogenases in H 2 photoproduction systems.

Overall, the literature is not fully consistent and difficult to assess because of the lack of precise report of experimental design and control experiments, regarding e.g. heating and photoelectric effects, or the details about the light sources and power, and the nature of the transparent materials the containers and electrochemical cells are made of. However, clear differences in photosensitivity have been observed between homologous FeFe hydrogenases, the enzyme from D. desulfuricans being particularly fragile in this respect. It is remarkable that this enzyme is also more O 2 sensitive than e.g. clostridial FeFe hydrogenases 14 , and sensitive to inhibition by sulfide whereas CpI hydrogenase apparently is not 70 . These variations are puzzling, considering that all these enzymes embed the same H-cluster, whose environment is very conserved. These observations demonstrate that structural details must influence the resistance to photochemical stress of FeFe hydrogenases. Such long distance effects are difficult to investigate. Regarding oxidative stress, we recently demonstrated that residues that are remote from the active site may change the flexibility of the protein and the reactivity in the 1st coordination sphere of the distal Fe. 89

Figure 1 -

 1 Figure 1 -Structure of the FeFe hydrogenase from C. pasteurianum CpI (white, PDB code: 6N59), C. reinhardtii HydA1 (red, PDB code: 3LX4 ) and D. vulgaris str. Hildenborough DdHydAB (yellow, PDB code: 1HFE). The FeS clusters and the H-cluster are shown in ball and stick representation.

  Katz et al. also recently observed by FTIR the photo-dissociation of exogenous CO from the H-cluster of CrHydA1 after blue light irradiation at 460 nm and 80 K, mimicking the resonance raman spectroscopy conditions, under which only the H ox and H red states result104 .This is a post-peer-review, pre-copyedit version of an article published inSustainable Energy & Fuels (2021) 10.1039/D1SE00756D Laser irradiation at low temperature of an H 2 -treated C. reinhardtii HydA1 sample induces an intramolecular electron transfer from the dinuclear cluster to the [4Fe4S] H subsite 104 . Lorent et al. recently studied the effect of light on the production of hydride-bound states under cryogenic conditions, in absence of CO 105 . The authors use photoactivation in combination with cryogenic infrared (IR) and electron paramagnetic resonance (EPR) spectroscopy to study the [FeFe]-hydrogenase from C. reinhardtii (CrHydA1), discovering two new terminal hydride species.

FIgure 3 -

 3 FIgure 3 -Species formed by irradiation of FeFe hydrogenases inhibited by CO at cryogenic temperatures, according to Chen67 and Roseboom66 . Irradiation of H ox -CO state at temperatures between 30 and 150 K produces essentially H ox , while no CO photolysis is observed at T>150 K due to rapid recombination of the CO ligand. Additionally, all light-induced changes can be completely reversed by warming the sample to 150 K in the dark.

  CO ligands. A difference is observed between C. pasteurianum, D. desulfuricans and C. reinhardtii hydrogenases: red light illumination exchanges the bridging CO of CpI and DdHydAB hydrogenases, whereas blue light is required to obtain the same effect in the case of CrHydA1. Stripp et al. propose the structures of the different forms from the IR spectra using a DFT model of the H-cluster with truncated cysteines, and larger models that include amino acids around the active site. Regarding the H ox -CO species, these calculations support the presence of the distal CN - This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D

Figure 5 .

 5 Figure 5. General mechanisms of catalysis and photocatalysis by [FeFe]-hydrogenase binuclear models.

  spectra of D. desulfuricans HydAB initially in the H ox state show This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy & Fuels (2021) 10.1039/D1SE00756D

Table 1 .

 1 Summary of the literature about the effect of light on distinct FeFe hydrogenases, under various conditions of irradiation: wavelength λ, power P, temperature (room temperature, 'RT', or cryogenic conditions, 'Cryo').

	Sustainable Energy

  According to time-dependent DFT (TDDFT) calculations, the H ox UV-visible spectrum is characterized by the presence of charge-transfer excited states

	This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy
	& Fuels (2021) 10.1039/D1SE00756D

  Besides the design and investigation of photocatalytic systems, an interest has grown in photochemistry of the simple Fe 2 S 2 model. Starting with the investigation of Fe 2 (CO) 6 (pdt) in 2008,133 many models with the general formula Fe 2 (CO) 4 (L) 2 (bridging ligand) have been studied (Table2). The synthetic Fe(CO) 4 (CN) 2 (adt) compound (used for the in vitro reconstitution of [FeFe]-hydrogenases) is catalytically inert155,156 , but the closely related

	This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy This is a post-peer-review, pre-copyedit version of an article published in Sustainable Energy
	& Fuels (2021) 10.1039/D1SE00756D & Fuels (2021) 10.1039/D1SE00756D
			Fe 2 (L) 6 (bridging ligand)
	Reference	λ(nm)	Ligands	bridging ligand
	Ridley et al., 2008 133 Fe(CO) 4 (CN) 2 (pdt) is active in H 2 photoproduction using [Ru(bpy) 3 ] 2+ as photosensitizer and 350 CO pdt
	Brown-McDonald et al., 2009 134 sodium ascorbate as sacrificial electron donor, upon irradiation with visible light (λ > 400 365 CO pdt
	Bertini et al., 2009 135 nM) 156,157	-	CO	pdt
	Strewart et al., 2010 136	350	CO	pdt
	Kaziannis et al., 2010 137	355,435	CO	pdt
	Marhenke et al., 2011 138	365	CO	pdt
	Stromberg et al., 2011 139	-	CO	pdt,edt
	Bingaman et al., 2012 140	266,355, 532	CO	pdt
	Muhammad et al., 2012 141	355	CO	pdt
	Kania et al., 2012 142	355	4CO, 2PMe 3	pdt
	Johnson et al., 2013 143	355,532	4CO, 2PMe 3	pdt
	Bertini et al., 2014 144	-	CO	pdt
	Caplins et al., 2014 145	400	CO	bdt
	Thornley et al., 2015 146	330, 400	CO	pdt
	Goy et al., 2016 147	-	CO	fluorene; silafluorene
	Hunt et al., 2016 148	365	4CO,2CN -	pdt
	Bertini et al., 2018 149	-	4CO, 2PMe 3	edt
	Meyer et al., 2018 150	266,355, 532	4CO, 2PMe 3	edt
	Stromberg et al., 2018 151	266, 400	4CO, 2CN -	pdt
	Thornley et al., 2019 152	200, 525	4CO, bipy/phen	pdt
	Bertini et al., 2021 153	-	CO, μCO	pdt
			mono substituted
	Meyers et al., 2021 154	400	CN -	edt,pdt
	Table 2 -Summary of the literature about the photochemistry of the Fe 2 S 2 model complexes
	(pdt=propandithiolate ; edt=etandithiolate ; bipy=2,2′-bipyridine ; phen=1,10-phenanthroline ; PMe 3 =
	trimethylphosphine).

Acknowledgements

The work in Marseille was supported by CNRS, Aix Marseille Université and the Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme (ANR-11-IDEX-0001-02). The French authors are part of the French bioinorganic chemistry network (www.frenchbic.cnrs.fr).


Xuewen Qian

Fellow, IEEE Marco Di Renzo
email: marco.di-renzo@universite-paris-saclay.fr

M X Di Qian

Renzo

K-Means Clustering-Aided Non-Coherent Detection for Molecular Communications

Keywords: Molecular communications, inter-symbol interference, channel state information, non-coherent detection, K-means clustering. I

In this paper, we consider non-coherent detection for molecular communication systems in the presence of inter-symbol-interference. In particular, we study non-coherent detectors based on memorybits-based thresholds in order to achieve low bit-error-ratio (BER) transmission. The main challenge of realizing detectors based on memory-bits-based thresholds is to obtain the channel state information based only on the received signals. We tackle this issue by reformulating the thresholds through intermediate variables, which can be obtained by clustering multi-dimensional data from the received signals, and by using the K-means clustering algorithm. In addition to estimating the thresholds, we show that the transmitted bits can be retrieved from the clustered data. To reduce clustering errors, we propose iterative clustering methods from one-dimensional to multi-dimensional data, which are shown to reduce the BER. Simulation results are presented to verify the effectiveness of the proposed methods.

the receiver [START_REF] Farsad | A novel experimental platform for in-vessel multi-chemical molecular communications[END_REF] and [START_REF] Qian | Receiver design in molecular communications: An approach based on artificial neural networks[END_REF]. However, ML-based receivers must be trained with long sequences of symbols. The performance of ML-based receivers, in addition, highly depends on the quality of the data used for training [START_REF] Zappone | Model-aided wireless artificial intelligence: Embedding expert knowledge in deep neural networks for wireless system optimization[END_REF].

B. Related Works

Non-coherent detection methods for MC systems are discussed in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF], where the authors obtain closed-form expressions of the optimal detection metric for a multi-symbol detector that uses an approximated probability mass function for the Poisson channel. In addition, the authors propose a blind single-symbol detector based on a constant threshold. However, the authors assume no ISI in the system model which may hold true when the symbol length is long enough and the symbol rate is low. In [START_REF] Li | Local convexity inspired low-complexity noncoherent signal detector for nanoscale molecular communications[END_REF], the authors consider the impact of ISI and exploit the local convexity of the received signals in order to detect the symbols from the difference of molecular concentration. To capitalize on the local convexity of the received signals, the transmitter needs, however, to release a sufficiently large number of molecules. In [START_REF] Li | Csi-independent non-linear signal detection in molecular communications[END_REF], a non-linear detector with an adaptive threshold is analyzed. The approach is based on the quick-rising and slow-decaying trends of the received signals after passing through a filter. In [START_REF] Liu | Unsupervised clustering-based non-coherent detection for molecular communications[END_REF], the authors exploit the unsupervised fuzzy C-mean approach to detect symbols from the quick-rising and slow-decaying of the processed signals.

In [START_REF] Luo | Non-coherent signal detection technique for mobile molecular communication at high data rates[END_REF], the authors propose an approach based on the energy difference of the received signals.

This method requires large numbers of molecules in order to exploit the energy difference.

Non-coherent detection is investigated in Poisson channels for application to optical wireless communications in [START_REF] Gong | Channel estimation and signal detection for optical wireless scattering communication with intersymbol interference[END_REF]. However, current non-coherent detection schemes assume negligible ISI. In summary, the aforementioned approaches either do not consider ISI, or require large numbers of molecules and, in general, cannot achieve low BER performance.

C. Contributions

In this paper, motivated by the results in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF] and [START_REF] Qian | Molecular communications: Model-based and data-driven receiver design and optimization[END_REF], we devise multi-memory-bit thresholdbased schemes without prior CSI for the reliable transmission of data in the presence of ISI.

The specific novelty and contributions made by this paper can be summarized as follows.

• We reformulate the multi-memory-bit threshold-based detectors via intermediate variables, i.e., the average number of received particles as a function of the memory bits. Unlike threshold-based methods that use the prior knowledge of the CSI, the considered intermediate variables can be obtained by applying clustering methods to the received signals.

• In order to reduce the clustering error, we construct multi-dimensional data from the received signals and designate the initial centroids of the clusters using the largest received signal. To the best of the authors' knowledge, it is the first time in the communication literature that the main features are extracted by increasing the dimension for clustering data. Conventional methods, on the other hand, are based on decreasing the dimension for clustering data [START_REF] Xu | Survey of clustering algorithms[END_REF]. In addition, we can infer the transmitted symbols from the clustered data. We design and analyze two approaches, which are referred to as direct clustering-based inference and clustering-plus-threshold detection.

• We devise an iterative method to further improve the BER performance and to avoid clustering failures in high data rate systems, i.e., in the presence of non-negligible ISI.

The essence of the proposed method is to iteratively apply clustering from one-dimensional data to high-dimensional data. In each iteration, we construct the initial centroids based on the estimated centroids in the previous iteration. By using the proposed iterative method, we show that the clustering errors can be reduced.

D. Paper Organization

The remainder of this paper is organized as follows. In Section II, the system model is introduced. In Section III, the multi-memory-bit threshold is reformulated via intermediate variables. In Section IV, the proposed non-coherent detection method based on multi-dimensional clustering is introduced. In Section V, an iterative method for estimating the initial centroids for application to multi-dimensional clustering is introduced. In Section VI, numerical results and the computational complexity of the proposed schemes are discussed. Finally, Section VII concludes this paper.

II. SYSTEM MODEL

We consider a three-dimensional unbounded MC system without flow that consists of a point transmitter and a spherical absorbing receiver, as depicted in Fig. 1. We assume that each information particle diffuses randomly and independently through the medium. The particles are assumed not to degrade rapidly, thus resulting in ISI. We assume that the temperature is constant and the viscosity remains unchanged during the whole transmission. Thus, the diffusion coefficient D [START_REF] Hiyama | Molecular Communication[END_REF] is assumed to be constant. Assuming that the transmitter is located in a = (0, 0, 0) and the receiver is located in b = (b x , b y , b z ), the hitting rate is expressed as follows [START_REF] Damrath | Low-complexity adaptive threshold detection for molecular communication[END_REF], [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF]:

f hit (t) = r(d -r) d √ 4πDt 3 e -(d-r) 2 4Dt (1)
where ab = d is the distance between the center of the transmitter and the center of the receiver, and r denotes the radius of the receiver. A binary CSK modulation scheme is assumed1 , and the bit transmitted in the ith time slot is denoted by s i . During the ith slot, the transmitter releases N T X information particles if s i = 1, otherwise the transmitter does not release any particles. We assume that the N T X information particles are released in a very short time such that we can ignore the release time effect on the received signal. The hitting probability of an absorbing receiver in the (i -1)th time slot of duration T is:

P i-1 = iT (i-1)T f hit (t)dt = r d erfc( d -r √ 4DiT ) -erfc( d -r 4D(i -1)T ) (2)
where erf(y) = y 0 2 √ π e -x2 dx and erfc(y) = 1 -erf(y). Let C j = N T X P j denote the average number of received particles in the jth time-slot after the release of N T X particles. For clarity, C 0 denotes the average number of received particles in the symbol of interest. Thus, the number of received particles [START_REF] Mosayebi | Receivers for diffusion-based molecular communication: Exploiting memory and sampling rate[END_REF] in the ith time-slot follows a Poisson distribution 2 :

r i ∼ Poisson(I i + s i C 0 ) (3)
where

I i = λ0 T + L j=1
s i-j C j is the sum of the ISI and background noise, λ0 is the background noise power per unit time, and L denotes the length (memory) of the Poisson channel.

We define the signal-noise-ratio (SNR) as follows:

SNR = 10 log 10 C 0 2 λ0 T (4)
since the information bits are assumed to be equiprobable. Thus,

N T X = 2 λ0 T 10 SNR 10 P 0 (5)
Therefore, the probability of receiving r i information particles is:

Pr(r i |I i + s i C 0 ) = e -(I i +s i C 0 ) (I i + s i C 0 ) r i r i ! (6)

III. MULTI-MEMORY-BIT THRESHOLD REFORMULATION

To detect the symbols, existing methods are based on calculating thresholds based on previously detected bits, which are referred to as memory bits. The main idea is to exploit previously detected bits, even erroneously estimated, in order to detect new transmitted bits based on the number of received particles. If the number of particles is below the threshold, a binary zero is estimated. Otherwise, a binary one is estimated. We first discuss the case of knowing the exact channel length L, and then we consider the case of limited channel information. The analytical threshold based on the knowledge of L memory bits can be obtained from the following equation:

Pr(r i |s i = 0, s i-j , 1 ≤ j ≤ L) = Pr(r i |s i = 1, s i-j , 1 ≤ j ≤ L) (7)
where Pr(r i |s i = 0, s i-j , 1 ≤ j ≤ L) and Pr(r i |s i = 1, s i-j , 1 ≤ j ≤ L) denote the probabilities of receiving r i information particles in the ith time slot conditioned upon the previously transmitted symbols s i-j for 1 ≤ j ≤ L and the current symbol s i = 0 or s i = 1, respectively. These probabilities correspond to [START_REF] Jamali | Channel estimation for diffusive molecular communications[END_REF].

By solving [START_REF] Fang | Symbol-by-symbol maximum likelihood detection for cooperative molecular communication[END_REF] based on (6), the optimal threshold τ | s i-j ,1≤j≤L conditioned on the previously transmitted symbols s i-j for 1 ≤ j ≤ L is as follows:

τ | s i-j ,1≤j≤L = C 0 ln(1 + C 0 I i ) = C 0 ln(1 + C 0 λ0 T + L j=1 s i-j C j ) (8)
Since the symbols s i-j for 0 ≤ j ≤ L are unknown, the previously estimated symbols ŝi-j for 0 ≤ j ≤ L are used to compute the threshold in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF], i.e., τ

| ŝi-j ,1≤j≤L = C 0 ln(1+ C 0 λ0 T + L j=1 ŝi-j C j ) .
Based on the computed threshold, ŝi is demodulated as follows:

ŝi =    1 r i > τ | ŝi-j ,1≤j≤L 0 r i ≤ τ | ŝi-j ,1≤j≤L (9)
Besides necessitating the bits detected in the previous time slots, the computation of the threshold in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF] relies upon the prior knowledge of the CSI, i.e., the coefficients C j are assumed to be known. To avoid using prior information about the channel, i.e., C j for 0 ≤ j ≤ L, when calculating the threshold in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF], we resort to intermediate variables, i.e., the average number of received particles ri | s i-j ,0≤j≤L that is defined as follows:

ri | s i-j ,0≤j≤L = I i + s i C 0 = λ0 T + L j=0 s i-j C j (10)
Equation ( 10) corresponds to the theoretical average number of particles, and, therefore, it depends on the variables C j . In the following, we first show how to formulate the thresholds in [START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF] by exploiting the intermediate variables in [START_REF] Mosayebi | Type-based sign modulation and its application for isi mitigation in molecular communication[END_REF], instead of using C j explicitly. In this step the CSI is still required. Then, we elaborate how to estimate the intermediate variables without necessitating CSI.

A. Computing the Thresholds Using ri | s i-j ,0≤j≤L

In this subsection, we show how the thresholds can be formulated only in terms of average number of particles ri | s i-j ,0≤j≤L . In particular, in order to formulate the problem analytically and understand the rationale of the proposed approach, we consider, just for this section, that ri | s i-j ,0≤j≤L is given by its analytical expression in [START_REF] Mosayebi | Type-based sign modulation and its application for isi mitigation in molecular communication[END_REF].

The threshold in (8) can be rewritten, in an equivalent form, only in terms of ri | s i-j ,0≤j≤L , as follows:

τ | s i-j ,1≤j≤L = C 0 ln(1 + C 0 λ0 T + L j=1 s i-j C j ) = C 0 + λ0 T + L j=1 s i-j C j -( λ0 T + L j=1 s i-j C j ) ln( λ0 T +C 0 + L j=1 s i-j C j λ0 T + L j=1 s i-j C j ) = ri | s i =1,s i-j ,1≤j≤L -ri | s i =0,s i-j ,1≤j≤L ln( ri | s i =1,s i-j ,1≤j≤L ri | s i =0,s i-j ,1≤j≤L ) (11)
In practice, the memory length used by the receiver may be limited and may be smaller than L. By assuming that L ≤ L memory bits are used, the threshold can be obtained by solving the equation Pr(r i |s i = 0, s i-j , 1 ≤ j ≤ L) = Pr(r i |s i = 1, s i-j , 1 ≤ j ≤ L) where:

Pr(r i |s i-j , 0 ≤ j ≤ L) = 1 2 L-L s i-j ,L+1≤j≤L Pr(r i |s i-j , 0 ≤ j ≤ L) (12)
Due to the analytical complexity of (12) it is, however, difficult to compute an analytical expression for the detection threshold by imposing Pr(r i |s i = 0, s i-j , 1 ≤ j ≤ L) = Pr(r i |s i = 1, s i-j , 1 ≤ j ≤ L). To overcome this analytical complexity, we use an approximation for Pr(r i |s i-j , 0 ≤ j ≤ L). In particular, we approximate [START_REF] Kilinc | Receiver design for molecular communication[END_REF] with Pr approx (r i |s i-j , 0 ≤ j ≤ L) as follows:

Pr approx (r i |s i-j , 0 ≤ j ≤ L) = e -r i | s i-j ,0≤j≤L (r i | s i-j ,0≤j≤L ) r i r i ! (13)
which is obtained by averaging

I i = λ0 T + L j=1
s i-j C j in (6) with respect to s i-j for L+1 ≤ j ≤ L, which yields:

ri | s i-j ,0≤j≤L = L j=0 s i-j C j + L j=L+1 C j /2 + λ0 T (14)
By using Pr approx (r i |s i-j , 0 ≤ j ≤ L), the detection threshold can be computed by solving the equation Pr approx (r

i |s i = 0, s i-j , 1 ≤ j ≤ L) = Pr approx (r i |s i = 1, s i-j , 1 ≤ j ≤ L), which yields: τ | s i-j ,1≤j≤L = ri | s i =1,s i-j ,1≤j≤L -ri | s i =0,s i-j ,1≤j≤L ln( ri | s i =1,s i-j ,1≤j≤L ri | s i =0,s i-j ,1≤j≤L ) (15)
As anticipated, the threshold in ( 15) is a theoretical formulation for the time being. In practice, it needs to be obtained (estimated) only from the received data without any prior information. To this end, we replace ri | s i-j ,0≤j≤L with its estimate ri | s i-j ,0≤j≤L , which is introduced and defined in the next sub-section. Therefore, the detection threshold is expressed as follows:

τ | s i-j ,1≤j≤L = ri | s i =1,s i-j ,1≤j≤L -ri | s i =0,s i-j ,1≤j≤L ln( ri | s i =1,s i-j ,1≤j≤L ri | s i =0,s i-j ,1≤j≤L ) (16)
Since s i-j for 0 ≤ j ≤ L are unknown, the estimates ŝi-j are used instead, which results in the following formulation for the threshold:

τ | ŝi-j ,1≤j≤L = ri | ŝi =1,ŝ i-j ,1≤j≤L -ri | ŝi =0,ŝ i-j ,1≤j≤L ln( ri | ŝi =1,ŝ i-j ,1≤j≤L ri | ŝi =0,ŝ i-j ,1≤j≤L ) (17)
Based on this reformulation, each symbol is demodulated as follows:

ŝi =    0, r i ≤ τ | ŝi-j ,1≤j≤L 1, r i > τ | ŝi-j ,1≤j≤L (18)
B. Estimation of ri | s i-j ,0≤j≤L from the Received Data

As anticipated, our objective is to develop non-coherent detection schemes that do not need CSI to operate. For example, we cannot rely on the knowledge of C j for j = 0, ..., L to compute the detection thresholds. In this sub-section, we show an example of how ri | s i-j ,0≤j≤L can be obtained by only using the received data, which implies that the threshold is ( 17) can be obtained without prior CSI information.

As a case study, we show how the intermediate variables ri | ŝi-j ,0≤j≤L are obtained for L = 1.

From ( 12), the exact probability mass function of receiving r i particles given s i and s i-1 is:

Pr(r i |s i , s i-1 ) = 1 2 L-1 s i-j ,2≤j≤L e -( λ0 T + L j=0 s i-j C j ) ( λ0 T + L j=0 s i-j C j ) r i r i ! (19)
Theoretically, ri | s i ,s i-1 is defined and can be computed as follows:

ri | s i ,s i-1 = E[r i | s i ,s i-1 ] = ∞ r i =0 r i Pr(r i |s i , s i-1 ) (20)
Let us consider a sequence of received particles {r 1 , ..., r n , ..., r K } of length K. The theoretical average number of received particles in (20) can be estimated empirically from the sequence of K observations. In practice (empirically), the probability Pr(r i |s i , s i-1 ) can be interpreted as a ratio [32, Section 19.2]:

Pr(r i |s i , s i-1 ) = lim Ns i ,s i-1 →∞ N r i ,s i ,s i-1 N s i ,s i-1 (21)
where N r i ,s i ,s i-1 denotes the number of elements within the sequence of K observations that are equal to r i and for which

[s n , s n-1 ] = [s i , s i-1 ],
and

N s i ,s i-1 = ∞ r i =0 N r i ,s i ,s i-1 denotes the number of elements within the sequence of K observations for which [s n , s n-1 ] = [s i , s i-1 ]
(but the number of particles is not necessarily equal to r i ). Therefore, (20) can be rewritten as follows:

ri | s i ,s i-1 = lim Ns i ,s i-1 →∞ ∞ r i =0 r i N r i ,s i ,s i-1 N s i ,s i-1 = lim K→∞ K n=0 r n κ n,s i ,s i-1 K n=0 κ n,s i ,s i-1 (22)
where

κ n,s i ,s i-1 = 1 if [s n , s n-1 ] = [s i , s i-1 ] and κ n,s i ,s i-1 = 0 otherwise. The variable κ n,s i ,s i-1
can be regarded as an indicator variable that provides information on whether the observation r n belongs to a group of observations for which

[s n , s n-1 ] = [s i , s i-1 ].
In the literature [START_REF] Jain | Data clustering: a review[END_REF], a group of observations is referred to as a cluster, a single observation of the cluster is referred to as a point of the cluster, and the condition that defines the cluster, i.e., the condition

[s n , s n-1 ] = [s i , s i-1 ]
, is referred to as the label of the cluster. In addition, the arithmetic mean of all the points in a subset X of R l is referred to as the centroid of the cluster [START_REF] Bovik | Handbook of image and video processing[END_REF] and is defined, in general terms, as follows:

µ = xg(x)dx g(x)dx ( 23
)
where x is a vector in R l , the integrals are computed over the whole space R l , and g(x) is a characteristic function that can be defined for different purposes. Based on ( 22) and ( 23), ri | s i ,s i-1 can be viewed as the centroid of the cluster of observations whose label is

[s i , s i-1 ].
Thus ri | s i ,s i-1 is obtained by computing the centroid of the cluster with label [s i , s i-1 ] that is obtained from the set of received signals (the observations) {r 1 , ..., r n , ..., r K }.

Based on [START_REF] Sabu | Detection probability in a molecular communication via diffusion system with multiple fully-absorbing receivers[END_REF], we evince that the thresholds of interest can be computed by calculating the centroids in [START_REF] Farsad | A novel experimental platform for in-vessel multi-chemical molecular communications[END_REF]. However, [START_REF] Sabu | Detection probability in a molecular communication via diffusion system with multiple fully-absorbing receivers[END_REF] necessitates the knowledge of the symbols s n , which are usually unknown. Our objective is to estimate ri | s i ,s i-1 from the received signals without knowing the symbols s n . Estimating the centroids in [START_REF] Sabu | Detection probability in a molecular communication via diffusion system with multiple fully-absorbing receivers[END_REF] or [START_REF] Farsad | A novel experimental platform for in-vessel multi-chemical molecular communications[END_REF] without any prior information can be obtained by applying clustering methods, which partition the received signals into several groups, each one corresponding to a given label [s i , s i-1 ], and then obtaining the corresponding centroids.

In summary, the thresholds for data detection can be obtained by applying clustering methods directly to the received signals without any prior information on the CSI and without knowledge of the transmitted symbols. In particular, the thresholds can be retrieved from the centroids of the clusters. In the next section, we detail how to use clustering for non-coherent detection in MC systems, and, in particular, we introduce our approach based on multi-dimensional clustering.

IV. CLUSTERING-BASED NON-COHERENT DETECTION

With the aid of the reformulated thresholds in Section III, we show that data detection can be realized without prior CSI. In particular, we introduce non-coherent detection methods based on the K-means clustering algorithm. To this end, we first introduce some background information on clustering in general and the K-means algorithm in particular. To elucidate the operating principle of the proposed methods, we report some illustrations that are obtained by using the simulation setup in Table I (see Section VI) in the presence of mild ISI (T = 30∆T ) and severe

ISI (T = 20∆T

). It is worth mentioning that the algorithms proposed in this paper require the storage of symbols in order to be executed, similar to [START_REF] Mosayebi | Receivers for diffusion-based molecular communication: Exploiting memory and sampling rate[END_REF], [START_REF] Farsad | A novel experimental platform for in-vessel multi-chemical molecular communications[END_REF], [START_REF] Liu | Unsupervised clustering-based non-coherent detection for molecular communications[END_REF], [START_REF] Shlezinger | ViterbiNet: A deep learning based viterbi algorithm for symbol detection[END_REF]. At the time of writing, this may be difficult to be realized with biological circuits. The actual implementation of the proposed solutions by using biological circuits is an important issue that needs to be assessed in future research.

A. K-Means Clustering Algorithm

Assume that we have a data set {x 1 , x 2 , ...x N } of N observations and each element is a Ddimensional vector x n . The objective is to partition the data set into N c clusters whose centroids are D-dimensional vectors denoted by µ k . The centroid µ k is the mean of its clustered points. In clustering, κ n,k is a variable that represents if the distance between the nth observation and the kth centroid is smaller than the distance between that observation and the other centroids. The K-means algorithm is a clustering method that works iteratively in order to compute or estimate the centroids µ k and the indicator variables κ n,k from a set of observations. In particular, the K-means algorithm encompasses two steps that iteratively compute κ n,k and µ k at each step.

The K-means algorithm needs an initial estimate of the centroids to operate. The initial centroids µ k can be either randomly selected from the data set or other methods can be employed. The selection of the initial centroids is discussed in further text. The K-means clustering algorithm can be summarized as follows.

Step I: Assign x n to the closest cluster (initial, if this is the first iteration) centroid:

κ n,k =    1, if k = arg min j x n -µ j 2 0, otherwise (24)
Step II: Update the cluster centroid:

µ k = n κ n,k x n n κ n,k (25)
where n κ n,k is the number of points in the kth cluster.

The two steps are repeated until κ n,k does not change. In the next sub-section, we show how to apply this algorithm for non-coherent detection in MC systems.

B. Single-Dimensional Clustering: Challenges and Limitations

In ( 16), the empirical average number of received particles ri | s i-j ,0≤j≤L need to be as close as possible to the theoretical average number of received particles ri | s i-j ,0≤j≤L in order to obtain a low BER. In order to elucidate the complexity of the problem at hand, let us consider the example in [START_REF] Guo | Eavesdropper localization in random walk channels[END_REF] by assuming L = 1.

The objective is to estimate ri | s i-j ,0≤j≤L from the set of received observations, {r 1 , ..., r n , ..., r K }, that are obtained from (3) given the transmitted symbols {s 1 , ..., s n , ..., s K }. Since L = 1, we are interested in identifying four clusters in the received signals, which correspond to the labels [START_REF] Huang | Channel characterization for 1-D molecular communication with two absorbing receivers[END_REF] or equivalently [START_REF] Qian | Receiver design in molecular communications: An approach based on artificial neural networks[END_REF], we obtain:

[s i , s i-1 ] = [0, 0], [s i , s i-1 ] = [0, 1], [s i , s i-1 ] = [1, 0], and [s i , s i-1 ] = [1, 1]. From
ri | s i ,s i-1 = s i C 0 + s i-1 C 1 + L j=2 C j /2 + λ0 T (26)
The average number of received particles in ( 26) is the theoretical one that can be obtained from analysis. In order to understand the difficulty of obtaining ( 26) from the empirical data, we illustrate an example in Fig. 2. Based on a large set of empirical data (K = 2 16 samples), we calculate the empirical probability mass function (PMF) of the number of particles, by assuming that the labels [s i , s i-1 ] are known. How the labels can be estimated from the data is discussed in further text. From the empirical data, we apply the K-means clustering algorithm in [START_REF] Zappone | Model-aided wireless artificial intelligence: Embedding expert knowledge in deep neural networks for wireless system optimization[END_REF] and [START_REF] Li | Local convexity inspired low-complexity noncoherent signal detector for nanoscale molecular communications[END_REF] in order to estimate the centroids. The estimated average number of received particles ri | s i ,s i-1 is set equal to the estimated centroids. From Fig. 2, we observe that there is a nonnegligible difference between the theoretical values obtained from [START_REF] Li | Csi-independent non-linear signal detection in molecular communications[END_REF] and the empirical values estimated by using the K-means clustering algorithm, even if the labels are assumed to be known, i.e., the association between the centroids and the labels is error-free. The differences between ri | s i ,s i-1 and ri | s i ,s i-1 mainly originate from the overlapping areas of the conditional probabilities Pr(r i |s i , s i-1 ).

The simple example illustrated in Fig. 2 highlights the challenges of applying clustering for estimating the average number of received particles in MC systems and, therefore, the detection thresholds. In the following sub-section, we introduce our proposed approach for improving the estimation accuracy based on multi-dimensional clustering methods. 26)) and ri | si,si-1 corresponding to the centroids calculated from ( 24) and ( 25) and assuming that the labels [s i , s i-1 ] are known) for T = 30∆T

C. Multi-Dimensional Clustering

The proposed solution to increase the accuracy of estimating the average number of received particles from the empirical data relies on multi-dimensional clustering. The main idea is not to apply the K-means clustering algorithm in [START_REF] Zappone | Model-aided wireless artificial intelligence: Embedding expert knowledge in deep neural networks for wireless system optimization[END_REF] and [START_REF] Li | Local convexity inspired low-complexity noncoherent signal detector for nanoscale molecular communications[END_REF] to individual observations r n but to vectors of observations:

r n = [r n , r n-1 , ..., r n-L ] ( 27
)
where L is the number of memory bits that we can use for detection and for multi-dimensional clustering. In order to illustrate the advantages of the proposed approach, let us consider the same example as in Fig. 2, but by assuming that the K-means clustering algorithm in [START_REF] Zappone | Model-aided wireless artificial intelligence: Embedding expert knowledge in deep neural networks for wireless system optimization[END_REF] and ( 25) is applied to the vector [r n , r n-1 ], rather than to r n only.

The results are illustrated in Fig. 3, where the following notation is used. The cluster of points whose label is, e.g.,

[s i = 0, s i-1 = 1] is denoted by [r n , r n-1 ]| s i =0,s i-1 =1
. A similar notation is used for the other clusters. By definition, ri

| s i ,s i-1 = s i C 0 + s i-1 C 1 + L j=2 C j /2 + λ0 T and ri-1 | s i-1 = s i-1 C 0 + L j=1 C j /2 + λ0 T given the label [s i , s i-1 ]
. Therefore, the corresponding (a) Empirical distributions of [rn, rn-1]|s i ,s i-1 and theoretical centroids obtained from [START_REF] Luo | Non-coherent signal detection technique for mobile molecular communication at high data rates[END_REF]. In this case, the labels are assumed to be known a priori (b) Theoretical centroids obtained by using (28) (yellow circles) vs. estimated centroids obtained by using the Kmeans algorithms in ( 24) and ( 25) (blue circles). The labels are not known a priori but are obtained from the estimated centroids (the updated initial centroids with assigned labels) Fig. 3: Theoretical and estimated centroids obtained from clustering the data vector [r n , r n-1 ] (T = 30∆T ) theoretical centroids are computed from ( 14) as follows:

[r i , ri-1 ]| s i =0,s i-1 =0 = [ L j=2 C j /2, L j=1 C j /2] + [ λ0 T, λ0 T ] [r i , ri-1 ]| s i =0,s i-1 =1 = [ L j=2 C j /2, L j=1 C j /2] + [C 1 + λ0 T, C 0 + λ0 T ] [r i , ri-1 ]| s i =1,s i-1 =0 = [ L j=2 C j /2, L j=1 C j /2] + [C 0 + λ0 T, λ0 T ] [r i , ri-1 ]| s i =1,s i-1 =1 = [ L j=2 C j /2, L j=1 C j /2] + [C 0 + C 1 + λ0 T, C 0 + λ0 T ] (28)
Given a sequence of received signals (observations), {r 1 , ..., r n , ..., r K }, the vectors of obser-

vations are {[r 2 , r 1 ], ..., [r n , r n-1 ], ..., [r K , r K-1 ]}.
The empirical distributions of each cluster of points and the theoretical centroids (computed from [START_REF] Luo | Non-coherent signal detection technique for mobile molecular communication at high data rates[END_REF]) are depicted in Fig. 3(a). We apply the K-means clustering algorithm in [START_REF] Zappone | Model-aided wireless artificial intelligence: Embedding expert knowledge in deep neural networks for wireless system optimization[END_REF] and [START_REF] Li | Local convexity inspired low-complexity noncoherent signal detector for nanoscale molecular communications[END_REF] to the empirical vectors using some initial centroids µ k . The setup of the initial centroids is described in the next sub-section. In general, the exact labels [s n , s n-1 ] are unknown, which implies that the association between an estimated centroid and the corresponding label [s i , s i-1 ] is not known a priori. We solve this issue by appropriately selecting the initial centroids, and by deciding the association between centroids and labels at the beginning of the K-means algorithm. Based on our approach, the estimated final centroids inherit the labels associated to the initial centroids. If [s i , s i-1 ] is the label of the initial centroid µ k , then [s i , s i-1 ] will be the label of the estimated (final) centroid μk . Therefore, the choice of the initial centroids is important to ensure a correct labelling. By comparing Fig. 3(b) and Fig. 2, we observe that ri | s i ,s i-1 is closer to ri | s i ,s i-1 in Fig. 3(b), which results in better performance. Therefore, multi-dimensional clustering is shown to be more accurate than one-dimensional clustering.

D. Setup of the Initial Centroids

As mentioned, an important issue is to assign the correct labels to the initial centroids and to ensure that the association between the estimated centroids and the correct labels does not change when applying the K-means clustering algorithm. In order to motivate the proposed approach for the setup of the initial centroids and the corresponding labels, let us consider Fig. 4(a).

In the figure, the four black dots correspond to the points [0, 0], [0, max(r)], [max(r), 0], and

[max(r), max(r)], where max(r) returns that largest value of the observations {r 1 , ....r n ....r K }.

In our proposed approach, the initial centroids of the clusters are chosen to correspond to the four black dots in Fig. 4(a). The rationale for this choice lies in the fact that the maximum received signal can be readily computed from the empirical data. In addition, more importantly, due to the geometry of the problem at hand, our proposed choice leads to a low probability that the label of the centroids will change after the application of the clustering algorithm. By direct inspection of the figure, in fact, the four black dots provide us with the correct label of the centroids, i.e., [0, 0] corresponds to the label [s i , s i-1 ] = [0, 0], [0, max(r)] corresponds to the label [s i , s i-1 ] = [0, 1], etc. If K is sufficiently large, the proposed association between centroids and labels is expected not to change, with high probability, during the application of the K-means algorithm.

If the number of memory-bits L is greater than one, the initial centroids are assigned in a similar fashion. If L = 2, for example, the initial centroids are [0, 0, 0] whose label is

[s i , s i-1 , s i-2 ] = [0, 0, 0]; [0, 0, max(r)] whose label is [s i , s i-1 , s i-2 ] = [0, 0, 1]; [0, max(r), max(r)] whole label is [s i , s i-1 , s i-2 ] = [0, 1, 1]; etc.
Based on these proposed initial centroids and the corresponding labels, the K-means clustering algorithm is applied according to ( 24) and ( 25), which returns the final estimated centroids ( μk ) and the indicator variables κ n,k . In particular, κ n,k allows us to implicitly perform data detection at the end of the clustering process, since it informs us, by definition, if a vector of points belongs or not to a given cluster. Based on this remark, we propose two algorithms for data detection that are referred to as (1) direct clustering-based inference and (2) clustering-plus-threshold detection.

E. Direct Clustering-Based Inference and Clustering-plus-Threshold Detection

The proposed direct clustering-based inference algorithm is based only on clustering methods and does not exploit the thresholds in [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF] for data detection. The algorithm is given in Algorithm 1 and works as follows (assuming, e.g., L = 1). If the observation vector [r n , r n-1 ] belongs to the cluster with label [s i , s i-1 ], i.e., the corresponding indicators variables are κ n,k = 1, then the estimated bits are those of the corresponding label.

The proposed clustering-plus-threshold detection algorithm, on the other hand, combines together clustering methods and the estimated thresholds in [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF] for data detection. The algorithm is given in Algorithm 2 and it has one main difference compared with Algorithm 1: After finalizing the clustering process based on the K-means algorithm, the bits are not detected by using the indicator variables κ n,k , but the thresholds in [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF] are computed from the estimated centroids, which are then used for data detection by using [START_REF] Bao | Channel modeling of molecular communication via diffusion with multiple absorbing receivers[END_REF].

In order to appreciate the difference between the two algorithms and the advantages and limitations of clustering and using or not using the detection thresholds in [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF], we provide some results in Figs. 5678. These figures are obtained as follows. Consider, for example, the case study L = 1 for ease of exposition. After estimating the received data ŝn , all the pairs of observations Compute the detection thresholds using [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF] 7: Detect the symbols using [START_REF] Bao | Channel modeling of molecular communication via diffusion with multiple absorbing receivers[END_REF] [ŝ n , ŝn-1 ] = [START_REF] Nelson | Biological physics: Energy, information life[END_REF][START_REF] Nelson | Biological physics: Energy, information life[END_REF] are depicted in red, green, blue, and magenta colors, respectively. The estimated centroids are depicted as black dots. The figures report the empirical decision boundaries that are obtained after clustering (Algorithm 1) and after applying the empirical thresholds estimated from the centroids (Algorithm 2). We note that the two algorithms provide different results, and their advantages and limitations depend on the severity of the ISI and number of memory bits. It is worth mentioning that in Figs. 6 and8, due to the many memory bits used, there exist more than four clusters and the dimension of the clusters is greater than two. To be able to illustrate the results, we merge together the clusters whose L-tuples [r n , r n-1 , . . . , r n-L ] have the same pair

[ŝ n , ŝn-1 ]
. By using this approach, we obtain again four clusters that can be readily represented.

The merged clusters are visualized by using the same color-based code as for L = 1.

By direct inspection of the figures, the following conclusions can be drawn. From Fig. 5, i.e., for mild ISI, we observe that Algorithm 1 outperforms Algorithm 2 if L = 1. In particular, Fig. we evince that Algorithm 1 with few memory bits may be a sufficiently good approach for noncoherent data detection under mild ISI. Under severe ISI, on the other hand, the performance trends are different, as illustrated in Figs. 7 and8. We observe, in particular, that the best estimation performance is obtained by using Algorithm 2 with a large number of memory bits, The initial centroids obtained from max(r), e.g., [0, 0], [0, max(r)], [max(r), 0] and [max(r), max(r)] if L = 1, are, in general, not always sufficiently close to the theoretical centroids.

This may result in detection errors. In this section, therefore, we introduce a more accurate iterative-based approach for estimating the initial centroids.

In order to illustrate the proposed approach, let us consider L = 1 (two-dimensional clustering) and the cluster whose label is

[s i , s i-1 ] = [0, 0]. Based on (28), the theoretical centroid that corresponds to [s i , s i-1 ] = [0, 0] is [r i | s i =0,s i-1 =0 , ri-1 | s i =0,s i-1 =0
], where:

ri-1 | s i =0,s i-1 =0 = ri-1 | s i-1 =0 = L j=1 C j /2 + λ0 T (29) ri | s i =0,s i-1 =0 = L j=2 C j /2 + λ0 T (30)
From ( 29) and ( 30), we observe that a good estimate for the initial centroid is [a, a] with a = L j=1 C j /2 + λ0 T , since this estimate would be close to the theoretical centroid, and, in general, it is closer than the initial centroid [0, 0] that was used in the previous section. Even though a depends on C j for j = 1, 2, . . . , L, an estimate for it can be obtained from onedimensional clustering without any prior CSI. In particular, a is, by definition, approximately

equal to a = ri | s i =0 ≃ ri | s i =0
, where ri | s i =0 is the estimated centroid obtained from onedimensional clustering. A similar approach can be used to estimate the centroids of the other clusters. The new initial centroids constructed by using ri | s i are depicted in Fig. 4(b) for completeness. By direct inspection, we can infer that they are closer to the theoretical centroids, and, therefore, are expected to provide better performance.

In general, the proposed approach lies in setting the initial centroids for L-dimensional clustering from the estimated centroids obtained by applying (L -1)-dimensional clustering.

In turn, the initial centroids for (L -1)-dimensional clustering are obtained from the estimated centroids obtained by aplying (L -2)-dimensional clustering. This procedure can be iterated until one-dimensional clustering, whose two initial centroids can be initialized to 0 and max(r).

The proposed approach for setting the initial centroids in an iterative-based manner is summarized in Algorithm 3 and Algorithm 4 for application to iterative clustering-based inference and iterative clustering-plus-threshold detection, respectively. In particular, the proposed approach for setting the initial centroids works in an iterative manner, as follows: (i) one-dimensional clustering is applied to the observations {r 1 , ..., r K } by setting 0 and max(r) as the initial centroids; (ii) the K-means clustering algorithm is applied and the estimated centroids ri | s i =0 and ri | s i =1 are obtained; (iii) the estimated centroids ri | s i =0 and ri | s i =1 are used to construct new initial centroids for application to two-dimensional clustering. In particular, the four initial

centroids are set to [ ri | s i =0 , ri | s i =0 ], [ ri | s i =0 , ri | s i =1 ], [ ri | s i =1 , ri | s i =0 ], and [ ri | s i =1 , ri | s i =1 ]; (iv)
two-dimensional clustering is applied to the vectors of observations {[r 2 , r 1 ], ..., [r K , r K-1 ]} by using the initial estimated centroids; (v) the K-means (two-dimensional) clustering algorithm is applied again and the estimated two-dimensional centroids [ ri , ri-1 ]| s i ,s i-1 are obtained; (vi) the procedure is iterated until the clustering dimension L + 1.

In general terms, the initial centroids for application to the (l + 1)-dimensional clustering can be constructed from the estimated centroids obtained from l-dimensional clustering. In mathematical terms, let us denote the estimated centroids from l-dimensional clustering by [ rm , ..., rm-l+1 ]| s m-j ,0≤j≤l-1 . The initial centroid µ k that corresponds to the label [s i , ..., s i-l+1 , s i-l ] for (l + 1)-dimensional clustering can be obtained as follows:

µ k = [ rm | s m-j =s i-j ,0≤j≤l-1 , rm | s m-j =s i-1-j ,0≤j≤l-1 , rm-1 | s m-j =s i-1-j ,1≤j≤l-1 , ..., rm-l+1 | s m-l+1 =s i-l ] (31)
which does not require any prior CSI information. Set the new initial centroids to µ k by using [START_REF] Yilmaz | Three-dimensional channel characteristics for molecular communications with an absorbing receiver[END_REF] 6: end for 7: Obtain ri | ŝi-j ,0≤j≤L from the estimated centroids μk 8: Compute the detection thresholds using [START_REF] Noel | Joint channel parameter estimation via diffusive molecular communication[END_REF] 9: Detect the symbols using [START_REF] Bao | Channel modeling of molecular communication via diffusion with multiple absorbing receivers[END_REF] By considering the same case study as in Fig. 8, we illustrate the performance obtained by employing the proposed iterative approach in Algorithms 3 and 4 in Fig. 9. We observe that, even in the presence of severe ISI, much better estimation performance is obtained by a more accurate initial estimate of the initial centroids. These results confirm the effectiveness of the proposed iterative-based approach.

VI. NUMERICAL RESULTS

In this section, we report some simulation results in order to quantitatively analyze the performance of the proposed non-coherent detectors. In particular, we compare the proposed methods with a few benchmark schemes available in the literature. The simulation parameters are listed in Table I.

A. Computational Complexity

Before analyzing the error performance of the proposed clustering-based algorithms, we study the computational complexity of the proposed schemes in Table II. In our study, we assume that

Step I and Step II of the K-means clustering algorithms are repeated P times. For instance, if we consider K received signals stored and L memory bits, (the centroid dimension is L + 1 and there are 2 L+1 initial centroids), in Step I of the K-means clustering algorithm, the evaluation of the distance between one point and a centroid requires 2L+1 additions and L+1 multiplications.

For all points, thus, we obtain 2 L+1 K(2L + 1) additions and 2 L+1 K(L + 1) multiplications. In

Step II, we need (K -2 L+1 )(L + 1) additions and 2 L+1 multiplications. Since the computational complexity of Step II is negligible compared with the computational complexity of Step I, we omit the complexity of Step II and report only the computational complexity of Setp I in Table II.

Algorithm 3, 4 P K L+1 i=1 2 i (2i -1) P K L+1 i=1 2 i (2i -1)

B. Channels in the Presence of ISI

In order to evaluate the efficiency of the proposed schemes, we compare the BER performance of the proposed algorithms with some benchmark schemes, including: (i) a receiver that first estimates the channel by using the method proposed in [START_REF] Jamali | Channel estimation for diffusive molecular communications[END_REF] by using K train pilot symbols and then detects the bits by using the threshold in [START_REF] Qian | Molecular communications: Model-based and data-driven receiver design and optimization[END_REF]; (ii) a receiver that detects the symbols by using the threshold in [START_REF] Qian | Molecular communications: Model-based and data-driven receiver design and optimization[END_REF] by assuming that the CSI is perfectly known; (iii) a sequence-based receiver [START_REF] Shlezinger | ViterbiNet: A deep learning based viterbi algorithm for symbol detection[END_REF] that detects the symbols by assuming perfect CSI (including the knowledge of L) and by using a multi-symbol maximum likelihood (ML) detector based on (13); (iv) the ViterbiNet sequence-based receiver in [START_REF] Shlezinger | ViterbiNet: A deep learning based viterbi algorithm for symbol detection[END_REF], which knows L and is trained using 2 10 pilot symbols. In order to show how the channel estimation accuracy affects the detection performance, we evaluate the benchmark scheme for K train = 2 6 and K train = 2 10 . In Fig. 10 (mild ISI) and Fig. 11 (severe ISI), we compare the performance of the proposed algorithms against the considered benchmark schemes by assuming a different number of memory bits L and the case studies in which the CSI is unknown (practical case) and known (benchmark case). We observe that the ML-based sequence detector with perfect CSI knowledge yields the best BER performance, as expected, and provides us with a lower bound. The proposed Algorithm 3 and Algorithm 4 yield good performance if the number of memory bits is not too small, depending on the severity of the ISI. By assuming that the CSI is known (e.g., the initial centroids are perfectly known in our algorithms), we observe (see Fig. 10(c) and Fig. 10(d)) that the proposed algorithms provide good performance. This shows the robustness of the proposed approach, which works well even under known CSI. We note, in addition, that the ViterbiNet shows an error-floor, in general, because it is derived by assuming a Gaussian white noise, which is different from the considered model.

In particular, we observe that Algorithm 1 and Algorithm 2 do not perform well, in general, in the presence of ISI because of errors in estimating the centroids. By comparing Fig. 10(c) and Fig. 10 (d), we evince that Algorithm 1 and Algorithm 2 have a worse BER if L increases. This is due to the high dimension of the clustered data. On the other hand, Algorithm 3 and Algorithm 4 compensate for this issue thanks to the iterative approach that is used to estimate the centroids. Based on the computational complexity that can be afforded, which depends on L, Algorithm 3 and Algorithm 4 are the most suitable options. In particular, Fig. The same trends are observed in Fig. 11. In the presence of severe ISI, it is necessary, however, to use a large number of memory bits L in order to outperform the ViterbiNet, as illustrated in Fig. 11(a).

C. Channels Without ISI

In order to further assess the robustness of the proposed algorithms, we evaluate the BER performance under ISI-free channels. In particular, we compare the BER performance of our algorithms against the benchmark in [START_REF] Jamali | Constant-composition codes for maximum likelihood detection without csi in diffusive molecular communications[END_REF], which is based on a non-coherent detector that uses constant-composition codes. In addition, we consider a benchmark that assumes that the CSI is known and uses the theoretical threshold in [START_REF] Qian | Molecular communications: Model-based and data-driven receiver design and optimization[END_REF]. The results are illustrated in Fig. 12. We observe that Algorithm 1 gives similar performance as Algorithm 3, and Algorithm 2 gives similar performance as Algorithm 4. This is because no ISI is present and the initial estimate of the centroids is almost correct. In particular, Algorithms 2 and 4 outperform Algorithms 1 and 3 because the decoding thresholds are exploited. In addition, Algorithm 2 and Algorithm 4 yield a BER that is close to the benchmark scheme with known CSI, and they outperform the benchmark in [START_REF] Jamali | Constant-composition codes for maximum likelihood detection without csi in diffusive molecular communications[END_REF].

D. Extension to Multi-Level CSK Modulation

To explore the potential applications of the proposed clustering-based methods, we generalize the proposed approach to multi-level CSK modulation [START_REF] Galmes | Performance analysis of diffusion-based molecular communications with memory[END_REF]. As an example, we consider a four- level CSK modulation scheme. By using similar analytical steps as those for the derivation of ( 15), e.g., for one memory symbol, we obtain the threshold:

τ k | s i-1 = ri | s i =k,s i-1 -ri | s i =k-1,s i-1 ln( ri | s i =k,s i-1 ri | s i =k-1,s i-1 ) (32)
by imposing Pr approx (r i |s i = k -1, s i-1 ) = Pr approx (r i |s i = k, s i-1 ). If r i ≤ τ 1 | s i-1 , then ŝi = 0;

if τ 1 | s i-1 < r i ≤ τ 2 | s i-1 , then ŝi = 1; if τ 2 | s i-1 < r i ≤ τ 3 | s i-1
, then ŝi = 2; otherwise, ŝi = 3.

In a 4-level CSK modulation scheme, there exist 16 clusters for the points [r n , r n-1 ] as depicted in Fig. 13(b). In this case, we consider Algorithms 3 and 4 for detection, and, in the first iteration, the initial centroids are set to µ 0 = 0, µ 1 = max(r)

3

, µ 2 = 2 max(r) 3 and µ 3 = max(r).

The numerical results are illustrated in Fig. 13. We observe that, in the considered case study,

Fig. 1 :

 1 Fig. 1: Three-dimensional MC system consisting of a point transmitter and a spherical absorbing receiver

Fig. 2 :

 2 Fig.2:The positions of ri | si,si-1 (calculated via[START_REF] Li | Csi-independent non-linear signal detection in molecular communications[END_REF]) and ri | si,si-1 corresponding to the centroids calculated from (24) and (25) and assuming that the labels [s i , s i-1 ] are known) for T = 30∆T

Fig. 4 :

 4 Fig. 4: The customized initial centroids and the corresponding clusters (T = 30∆T )

[r n , r n- 1 ] 1 Algorithm 2

 112 that are detected as [ŝ n , ŝn-1 ] = [0, 0], [ŝ n , ŝn-1 ] = [0, 1], [ŝ n , ŝn-1 ] = [1, 0], and Algorithm Direct clustering-based inference First Step: Clustering 1: Set L 2: Construct the data r n via (27) 3: Construct the initial centroids µ k 4: Cluster the r n using the K-means algorithm with the initial centroids µ k Second Step: Inference 5: Infer ŝn from the indicator variables κ n,k Clustering-plus-threshold detection First Step: Clustering 1: Set L 2: Construct the data r n via (27) 3: Construct the initial centroids µ k 4: Cluster the r n using the K-means algorithm with the initial centroids µ k Second Step: Detection 5: Obtain ri | ŝi-j ,0≤j≤L from the estimated centroids μk 6:

Fig. 5 :

 5 Fig. 5: Distribution of [r n , r n-1 ]| ŝn,ŝn-1 , estimated centroids and clustering boundaries/thresholds (L = 1 and T = 30 ∆T ).

Fig. 6 :

 6 Fig. 6: Distribution of [r n , r n-1 ]| ŝn,ŝn-1 and estimated centroids (L = 4 and T = 30 ∆T ).

Fig. 7 :

 7 Fig. 7: Distribution of [r n , r n-1 ]| ŝn,ŝn-1 , estimated centroids and clustering boundaries/thresholds (L = 1, T = 20 ∆T ).

Fig. 8 :

 8 Fig. 8: Distribution of [r n , r n-1 ]| ŝn,ŝn-1 and estimated centroids (L = 4, T = 20 ∆T ).

  i.e., L = 4 in Fig.8(b). Errors are still clearly visible, but the estimated clusters of points are closer to the expected ones, i.e., those shown in Fig.7(a). It is worth mentioning that in the presence of severe ISI the combination of clustering methods and detection thresholds provides better performance than using only clustering. The use of detection thresholds, in particular, allows us to correct the increased number of mis-detections that originate from increasing the dimension of clustering. These qualitative conclusions that are drawn from the direct inspection of the estimated clusters from the detected bits are corroborated in Section VI with the aid of BER simulations.V. ITERATIVE ALGORITHM FOR COMPUTING THE INITIAL CENTROIDSIn the previous section, we have illustrated the core ideas of the proposed clustering-based algorithms for non-coherent detection in MC systems, which may or may not use detection thresholds. The comparison of, e.g., Fig.7(a) and Fig.8(b) reveals, however, that detection errors still exist, especially in the presence of non-negligible ISI. The sources of error include the initial estimates of the centroids and the resulting estimates of the detection thresholds. In this section, motivated by these considerations, we introduce improved, iterative-based, clusteringbased methods in order to enhance the detection performance in the presence of severe ISI.

Algorithm 3 5 :Algorithm 4 4 :

 3544 Iterative clustering-based inference 1: Set the initial centroids µ 0 = 0 and µ 1 = max(r) 2: for l = 1 to L + 1 do 3: Construct the data r n = [r n , ..., r n-l+1 ] 4: Cluster r n using the K-means algorithm with the initial centroids µ k Set the new initial centroids to µ k by using (31) 6: end for 7: Infer ŝn from the indicator variables κ n,k Iterative clustering-plus-threshold detection 1: Set the initial centroids µ 0 = 0 and µ 1 = max(r) 2: for l = 1 to L + 1 do 3: Construct the data r n = [r n , ..., r n-l+1 ] Cluster r n using the K-means algorithm with the initial centroids µ k 5:

Fig. 9 :

 9 Fig. 9: Distribution of [r n , r n-1 ]| ŝn,ŝn-1 and estimated centroids (L = 4, T = 20∆T ).

(d) L = 4 Fig. 10 :

 410 Fig. 10: BER comparison of the proposed algorithms (mild ISI)

=2 6 )Fig. 11 :

 611 Fig. 11: BER comparison of the proposed algorithms (severe ISI)

Fig. 12 :

 12 Fig. 12: BER comparison of the proposed algorithms (no ISI, T = 100 ∆T and L = 1)

Algorithm 3 4 Fig. 13 :

 3413 schemes are inherently more vulnerable to wrong estimations of the decoding thresholds. Therefore, Algorithm 3 outperforms Algorithm 4 in the considered case study. This behavior can be understood by comparing Fig.13(b) with Fig.13(c) and Fig.13(d), in which we evince that Algorithm 3 yields a smaller classification error compared with Algorithm 4.

TABLE I :

 I Simulation parameters

	Parameter λ0	Value 100s -1
	Receiver radius r	45 nm
	Distance d	500 nm
	Diffusion coefficient D Discrete time length ∆T 9 us 4.265 × 10 -10 m 2 /s
	Channel length L	5
	K	2 12

TABLE II

 II

		: Computational complexity
	Algorithm	Additions	Multiplications
	Algorithm 1, 2 2 L+1 P K(2L + 1) 2 L+1 P K(L + 1)

The generalization to multilevel CSK modulation is discussed in Section VI-C.

The number of received particles ri can be modelled accurately by using a Binomial distribution. For ease of analysis, ri is approximated with a Poisson distribution. This is a good approximation if NT X is large[START_REF] Jamali | Non-coherent detection for diffusive molecular communication systems[END_REF].