Farad Khoyratee

Thomas Benneteau

Agnès Tixier-Mita

Sylvain Saïghi

Timothée Levi

BIO-HYBRID EXPERIMENTS USING TUNABLE REAL-TIME BIOMIMETIC NEURAL NETWORK

INTRODUCTION

Neuroscience has been the object of many studies and new fields these last decades leading to an era where technology and biology could be mixed in order to find advanced solutions for curing diseases. Neuroprosthesis require to take into consideration the biological activity of the cells and the connection between machines and biological cells. Understanding the mechanism underlying of affected cells also require to use a tunable and real time system such as a Field Programmable Gate Array (FPGA) to explore different conditions and parameters. Here, we propose to study neurological diseases by using digital hardware system reproducing the spike activity of Fast Spiking (FS), Regular Spiking (RS), Intrinsically Bursting (IB) and Low-Threshold Spiking (LTS) class of cortical neurons [START_REF] Pospischil | Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons[END_REF]. We use Hodgkin and Huxley (HH) [START_REF] Hodgkin | A quantitative description of membrane current and its applications to conduction and excitation in nerve[END_REF] which is the most bio-plausible model.

METHODS

Mathematical statements using HH paradigm are complex but it is possible to reduce them by using mathematical and numerical functions like COordinate Rotation DIgital Computer (CORDIC) to perform hyperbolic functions, Euler methods to solve derivative functions and non-restoring divisions. The more the equations are reduced, the more the speed and space are increased. Thus, implementing other biological parameters like axonal delay, synaptic noise is possible. Synapses and short-term plasticity have also been implemented and allow the chip to compute biomimetic locomotion with Central Patter Generator (CPG). The integrated system is fully pipelined and is composed of a main unit managing RAM, computation units and the Digital to Analog Converter (DAC) which is used to output the APs. Moreover, Real-time modification of biological parameters like reversal potential, ionic conductance or membrane capacitance is allowed. The results have been verified by comparing biological data in [START_REF] Hodgkin | A quantitative description of membrane current and its applications to conduction and excitation in nerve[END_REF] and electrical activities from the FPGA. Then, cross-correlation equations have been used to evaluate the proximity rate of the results.

RESULTS

Mathematical equations have been significantly simplified. Resources are reduced by 80% from original model. Regular and simplified equations of biomimetic neural network are implemented on a FPGA. Electrical activities of different classes of cortical neurons can be simulated in biological time scale. Mathematical tools show strong correlation between hardware simulations and biological recordings (see figure 1). This biomimetic system is ready for bio-hybrid experiments with 'in vitro' neuron culture.

CONCLUSION

A biomimetic chip has been designed using mathematical and numerical technics for bio-hybrid experiments. The simplification of the equations allow future work on the implementation of multi-compartmental model which takes account of different part of cells, characterization of disorders effect on a neural networks. This system can be connected to an in-vitro culture with a pattern using a microfluidic technic (see figure 2).

KEYWORDS:

Biomimetic, Neural Network, Hodgkin-Huxley, FPGA, Bio-hybrid Experiments


Stephany Mai Nishikawa

Soo Hyeon Kim

Yoshiho Ikeuchi

Teruo Fujii

Timothée Levi

Ionic Biomimetic

Microstimulation

la diffusion de documents scientifiques de niveau

recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Neurodegenerative diseases are incurable and debilitating conditions that influence cognitive and/or motor functions in millions of people worldwide. Neuroprosthesis are used today to support the quality of life but have yet to improve in their power consumption and biocompatibility issues. For the future realization of neuroprosthesis, understanding the neurophysiological behaviors and the investigations on the interaction of neuronal cell assemblies is therefore essential. Here, we propose a novel microfluidic system to investigate the response of the neurons directly stimulated by the potassium ions in vitro. Unlike the conventional way of stimulating neurons by pipetting or sending electrical stimulus, this biomimetic stimulation system enables precise and rapid control of ion concentration by utilizing microfluidic pneumatic valves, which mimics physiological stimulation in a body. By using the system, we demonstrated stimulations of the neuron cells with potassium ions and monitored neural activities with a calcium-specific fluorescence indicator.

CONCEPT OF THE DEVICE

The device was inspired by nervous system which is a network of neurons. Each neuron is connected through synapses, to a large number of other neurons, through dendrites and axons. The axon transmits electrochemical pulses (including KCl, NaCl) to communicate between neighboring neurons by transmitting electrical signals based on ionic flow between intracellular and extracellular environments. To simulate this physiological behavior, a device made with Polydimethylsiloxane (PDMS) has two chambers intra-and extracellular chambers connected by channels with high and low concentrations of potassium chloride actuated by pneumatic microvalves. An ion selective permeable membrane (Sigma-Aldrich, Nafion® perfluorinated resin solution 5wt.%) is used for and the concentration of potassium ions can be changed.

RESULTS

The device with integrated cell dish for placing living neurons has a thinness of less than 5mm allowing the observation through microscope during calcium imaging. Our device can generate a biomimetic stimulation accordingly to the sequenced air pressure. The opening of the pneumatic microvalve was sequentially repeated after three pulses of 0.1s, 0.25s and 0.5s with each separated with 0.2s of closure. The matured neurons were stained with calcium specific fluorescence to observe the neural activity. Fluorescence response was successfully observed only during the opening of pneumatic microvalve.

CONCLUSIONS

We successfully demonstrated stimulations of the neuron cells with potassium ions using microfluidic pneumatic microvalves operated by air pressure controller. The neural activities were visible according to the programmed air pressure applied, capable of generating a biomimetic behavior. We believe that our novel system can be utilized for the fundamental investigation of the response of the neurons stimulated by well-modulated ion stimulus to understand the neurophysiological behaviors.


Farad Khoyratee

Luo Zhongyue

Carole-Anne Volette

Thomas Benneteau

Romain Beaubois

Corentin Lange

Kazuyuki Aihara

Teruo Fujii

Yoshiho Ikeuchi

Timothée Levi

NEUROMIMETIC STIMULATION APPLIED ON CEREBRAL ORGANOIDS

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

INTRODUCTION

Millions of people are affected by neurological disorders. Brain-Machine Interfaces (BMIs) and neuroprosthesis have been object of extensive research and may represent a valid treatment for diverse neurological diseases. The realization of neuroprostheses requires technologies to interact with neuronal cell assemblies in the nervous system, and to drive them into a desired state or to produce a specific behavior. A hardware set-up that can be used to connect with biological system is a Spiking Neural Network (SNN) system which implements biologically realistic neural network models, spanning from the electrophysiological properties of one single neuron up to network plasticity rules. In this work, we present a digital hardware implementation of a biomimetic SNN reproducing a Central Pattern generator (CPG) [START_REF] Ambroise | Real-time biomimetic Central Pattern Generators in an FPGA for hybrid experiments[END_REF] for stimulating organoids [START_REF] Lancaster | Cerebral organoids model human brain development and microcephaly[END_REF], [START_REF] Kawada | Generation of a Motor Nerve Organoid with Human Stem Cell-Derived Neurons[END_REF] which are three-dimensional assemblies of neurons developed from human stem cells in vitro. By connecting organoids mimicking our nervous system with electronic SNNs, we intend to develop a model system for a brain-electronics interface.

METHODS

Most of rhythmic movements are programmed by central pattern-generating networks that comprise neural oscillators. This CPG consists of two neurons, connected to each other via inhibitory synapses including short term plasticity. We implement real-time biomimetic central pattern generators (CPG) into digital hardware (FPGA). This CPG mimics the Leech heartbeat neural network system [START_REF] Kawada | Generation of a Motor Nerve Organoid with Human Stem Cell-Derived Neurons[END_REF]. This system is composed of Izhikevich neurons, synaptic and plasticity systems. Our implementation needs few resources, consumes low power, and can be easily embedded. Our system has been validated by comparing to biological data. The CPGs provide bursts of spikes with frequency, variability and spike width in the ranges commonly observed in nature, and thus are optimal for biomimetic stimulations. Cerebral organoids were generated as described before [START_REF] Lancaster | Cerebral organoids model human brain development and microcephaly[END_REF]. After loaded with Fluo-4 AM dye, cerebral organoids were stimulated with CPGs and changes of fluorescent intensity were analyzed with an inverted microscope equipped with a sCMOS camera.

RESULTS

Three different CPGs with different periods are used (1,7s, 2,9s and 5s). Outputs of the CPG are triggered with the Calcium imaging camera recording. The stimulation applied to the isolated organoid is a 4V pulse of 800µs. Each spike of the CPG occurs one stimulation on the organoid. We can notice that the living organoid dynamically reproduced the pattern of oscillation applied by the CPG.

CONCLUSIONS

This biomimetic CPGs are working in real-time and can send biomimetic stimulations on living organoid. We validate the stimulation using Calcium imaging response. This bio-hybrid system will be used for investigating network dynamics in organoid and for controlling living cells to activate in a specific temporal pattern [START_REF] Aa | A model of a segmental oscillator in the leech heartbeat neuronal network[END_REF].


Anne-Claire Eiler

Junichi Sugita

Satoshi Ihida

Hiroshi Toshiyoshi

Katsuhito Fujiu

Timothée Levi

Agnès Tixier-Mita

Bio- Electrical Signal

Bioelectrical Signal Analysis of Mouse Cardiomyocyte Culture recorded on Thin-Film-Transistor Sensor Arrays

Keywords: Thin-Film-Transistor Arrays, Bioelectrical Signal Processing, Electrophysiology, Cardiomyocytes

  L'archive ouverte pluridisciplinaire

Introduction

Cardiomyocytes are primary muscle cells derived from heart tissue that can generate and conduct bioelectrical signals for cardiac contraction and blood flow. A * The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan problem occurring in the cellular bioelectrical network can range from minor to fatal inconvenience. [START_REF] Farwell | Electrical heart disease: Genetic and molecular basis of cardiac arrhythmias in normal structural hearts[END_REF] The cells can retain their physiological functions and thus provide a useful in vitro model to look at the beating rate, the duration and the shape of the field potential. In vitro research of the general behavior of cardiomyocytes can thus help understanding arrhythmias, long Q-T syndrome and cardiotoxicity, As a result, in vitro study of cardiomyocytes represent a valuable tool for drug discovery and cardiac research. However, the mechanism of the cardiac bio-signaling network is still poorly understood.

light of this problem, this paper proposes the analysis of bioelectrical signals of cardiomyocytes measured on an array of microelectrodes using a new Thin-Film-Transistor (TFT) sensor array. Measurements of the bioelectrical potentials of neurons were already demonstrated with TFT arrays. [START_REF] Shaik | Extracellular neural stimulation and recording with a Thin-Film-Transistor (TFT) array device[END_REF] The data flow generated by large arrays must be compressed to envision compact data acquisition systems. Hence, the electrical signals have been analyzed using a MATLAB program developed for bioelectrical processing of electrogenic cells. The recorded signals were filtered for the detection of spikes, and grouped into clusters according to their similar features. Through this analysis, the experiments demonstrated the possibility of obtaining accurate spike sorting and analysis from extracellular recordings on TFT arrays.

Measurement Method

Thin-Film-Transistor Arrays

The TFT arrays were used for extracellular recordings in vitro of cardiomyocytes. Here, the TFT technology wellknown for Liquid Crystal Display is used for biological applications. [START_REF] Tixier-Mita | Review on thinfilm transistor technology, its applications, and possible new applications to biological cells[END_REF] The standard type of TFT array comes in a pattern of 150 x 150 transparent microelectrodes and is mounted on a printed circuit board (PCB). Microelectrodes are composed of indium tin oxide (ITO) with a size of 100 x 100 µm 2 . The array of microelectrodes is controlled by an array of TFTs, which are used for switching ON/OFF their respective microelectrodes. The TFTs are controlled by means of gate and source/drain lines. The columns of the array control the gates of the TFTs, while the rows control the sources. When a 12V DC voltage is applied to one gate line, all the microelectrodes connected to that line are activated. Then one or more source lines wereconnected to a measurement system for sensing. The bioelectrical signals of cardiomyocytes were measured using a Multi-Channel Measurement System (MCS USB-ME32-FAI-System) and optical observations were performed with an inverted microscope. Fig. 1 describes the working principle of TFT arrays and the experimental setup.

Culture of Cardiomyocytes

Cardiomyocytes were dissociated into single-cell suspension from neonatal mice by combining mechanical dissociation with enzymatic degradation of the extracellular matrix, which maintains the structural integrity of tissues. The neonatal hearts were enzymatically digested using a neonatal heart dissociation kit for mouse from Miltenyi Biotec and a dissociator was used for the mechanical dissociation steps. After dissociation, the sample was filtered to remove any remaining larger particles from the single-cell suspension, and the lysis of red blood cells was performed.

Cardiomyocytes were finally cultured for 3 days on the TFT array devices without surface treatment.

Bioelectrical Analysis

Embedded signal processing is an essential step in the development of recording instrumentation. Here, spike sorting algorithm was used. This data processing technique consists in identifying the cells that contribute to the signal recorded by each microelectrode, their number and their spiking activities. [START_REF] Beche | Real-time adaptive discrimination threshold estimation for embedded neural signals detection[END_REF] The identified basic functions are ( 1

Bioelectrical Signal Recording

In this study, the extracellular electrical potentials of cardiomyocytes were first recorded on 28x28 microelectrodes. The measured noise level was approximately 50 µV. The bioelectrical activity of the culture of cardiomyocytes was confirmed by optical visualization of the cell contraction using the inverted microscope (Fig. 3.). Here, a line of 4 microelectrodes was selected and data were extracted for the ensuing processing of the bioelectrical signals.

Bioelectrical Signal Processing

Bioelectrical signal processing of the data from several measurements was performed. Each analysis provided valuable information about the bioelectrical signals of the cardiomyocytes. In this study, the data from a decrease in temperature has been analyzed.

Filtered Data

Raw data were first filtered to remove undesired signals according to their frequency. This low-pass filter passes signals with a frequency lower than the selected cutoff frequency of 200 Hz. The sampling rate of the recorded data was 10 kHz. Fig. 4. shows an example of data obtained before and after filtering.

Spike Detection and Alignment

The action potentials of cardiomyocytes are then distinguished from the noise according to a predefined threshold. The spikes are then aligned with respect to the maximum absolute value of the detected signal as shown in Fig. 5.

Spike Intervals

For each microelectrode, the spike intervals are classified. Histogramsdisplay the number of spikes according to their period (Fig. 6.).

Clustering

Spikes were finally divided into a number of clusters such that spikes of the same group are more similar than those in other groups. To reduce sample space dimension, Principal Component Analysis (PCA) was used. This statistical technique provides features that are directions in the high-dimensional space. In this study, each spike became a point in a 3-dimension space.

Discussion

The transparent TFT array device can provide a high spatial resolution of cell culture activity by surpassing currently available MEAs in terms of density of microelectrodes, and its cm-sized measurement surface area. With regard to the bioelectrical analysis, the spike sorting program allowed to extract useful information about the bioelectric signals generated by cardiomyocytes. A decrease of the beating rate with the temperature (from +37°C environment to room temperature) at around -6.4 beats per minute/°C on each microelectrode is confirmed by the program. A possible explanation would be that the temperature drop depresses the speed of ion exchange, as it increases the permeability of the membrane to ions. The program also identified a strong irregularity of the peak-to-peak voltage amplitude with the decrease in temperature. However, one main cluster of spikes is identified for each microelectrodes and share equivalent waveform of the mean spikes. This waveform is highly similar to the typical cardiac action potential in electrocardiogram with the typical T wave associated with ventricular repolarization. A spike raster plot also revealed the coincidence of spike occurrence for each microelectrode. This confirms the synchronicity of the bioelectrical conduction among the cell culture, which was also observed in parallel with the microscope.

Conclusion

In this paper, the bioelectrical activity of cardiomyocytes was successfully measured on TFT sensor arrays and analyzed by using spike sorting technique. This analysis confirms that TFT arrays can efficiently detect the bioelectrical signals generated by cardiomyocytes. Combining this technique with deeplearning algorithms 5 could allow the in vitro identification of abnormal cardiac cell conduction and aid for drug screening.

Fig. 1 .

 1 Fig. 1. (A) Working principle of the transparent TFT substrate, with a close-up view of the array of 100 μm square microelectrodes; (B) Experimental setup

Fig. 2 .

 2 Fig. 2. Workflow of the functions performed for spike sorting 3. Results

Fig. 3 .

 3 Fig. 3. Culture of cardiomyocytes on TFT array device, 3 days after cell seeding.

Fig. 4 .

 4 Fig. 4. (A) Raw data of the total measurement time (7 min) and (B) close-up view of a single spike on filtered data.

Fig. 5 .

 5 Fig. 5. (A) Filtered data and threshold lines in red for spike detection, and (B) spike alignment one 1 microelectrode.

Fig. 6 .

 6 Fig. 6. (A) Histogram of the spike intervals detected on 1 microelectrode when the temperature of the culture chamber is at +37°C and (B) at room temperature (RT). In this study, the spike intervals increased with decrease of temperature: from 190 ms (5.3 Hz) at +37°Ct to 270 ms (3.7 Hz) at RT.

Fig. 7 .

 7 Fig. 7. (A) Mean spikes of each 4 clusters identified after analysis of 4 microelectrodes, and (B) alignment of the mean spikes from each cluster.

© The 2020 International Conference on Artificial Life and Robotics (ICAROB2020), Feb. 20-23, B-Con Plaza, Beppu, Oita, Japan

© The 2020 International Conference on Artificial Life and Robotics (ICAROB2020),Feb. 20-23, B-Con Plaza, Beppu, Oita, Japan


Yang Xia
email: xiayang@sat.t.u-tokyo.ac.jp

Timothée Levi

Takashi Kohno

Real-Time

Real-time pattern recognition implementation on FPGA in multi-SNNs

Keywords: SNN, STDP, DSSN, FPGA, Ethernet

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Neuromorphic systems are designed by mimicking or being inspired by the nervous system, which aims to realize robust, autonomous, and power-efficient information. There are three common methods to realize the neuromorphic circuits, which are software [START_REF] Hines | NEURON: a tool for neuroscientists[END_REF] , analog hardware [START_REF] Natarajan | Hodgkin-Huxley Neuron and FPAA Dynamics[END_REF][START_REF] Grassia | Tunable neuromimetic integrated system for emulating cortical neuron models[END_REF] and digital hardware 4,[START_REF] Grassia | Digital hardware implementation of a stochastic two-dimensional neuron model[END_REF] . Software-based cannot achieve real-time processing. Compared to analog circuits, digital implementations generally consume higher power but are more scalable because they are less susceptible to noise and fabrication mismatch. Another advantages of FPGA devices are configurability, portability, and low-cost. Silicon neuronal network (SNN) is a neuromorphic circuit that reproduces the electrophysiological activities in the nervous system focusing on the spiking dynamics in the neuronal cells and their transmission via the synapses. Their application includes bio-inspired information processing such as pattern recognition [START_REF] Qiao | A reconfigurable on-line learning spiking neuromorphic processor comprising 256 neurons and 128K synapses[END_REF][START_REF] Merolla | A million spiking-neuron integrated circuit with a scalable communication network and interface[END_REF] and associative memory [START_REF] Li | An FPGA-Based Silicon Neuronal Network with Selectable Excitability Silicon Neurons[END_REF] as well as neuro-prosthetic devices 9 , 10 .

Digital SNNs are expected to achieve a very large-scale network comparable to the human brain in the future exploiting the scalability of the digital circuits. In these 5 years, several very large-scale SNNs with one million neurons were developed [START_REF] Davies | Loihi: a neuromorphic manycore processor with on-chip learning[END_REF] . Spike-timing-dependent plasticity (STDP) is a wellknown rule for updating the synaptic efficacy in SNNs, which uses only local information. Many biological experiments found evidence for STDP process in the synapses [START_REF] Markram | Regulation of synaptic efficacy by coincidence of postsynaptic aps and epsps[END_REF][START_REF] Jacob | Spike-timing-dependent synaptic depression in the in vivo barrel cortex of the rat[END_REF] . In this paper, we report an implementation of digital SNN with online STDP learning on FPGA. The model and implementation of our SNN are explained in the next section and section 3, respectively. Then results and conclusion follows.

Architecture of the network model

DSSN model

Neuronal models need to be chosen taking into account the balance between the reproducibility of neuronal activity and computational efficiency. Integrate-and-fire (I&F)-based models are able to be implemented by compact hardware, but they lack reproducibility of complex neuronal dynamics. Ionic-conductance models have high-ability of reproducing neuronal activities, but cost massive computational resources. The DSSN model is a qualitative neural model [START_REF] Kohno | Digital spiking silicon neuron: concept and behaviors in GJ-coupled network[END_REF] , which was designed for efficient implementation in digital circuits. The simplest version of DSSN model supports the Class I and II cells in Hodgkin's Classification [START_REF] Hodgkin | The local electric changes associated with repetitive action in a non-medullated axon[END_REF] .The differential equations of DSSN model are as follows.

𝑑𝑣 𝑑𝑡 = Ф 𝜏 (𝑓(𝑣) -𝑛 + 𝐼 0 + 𝐼 𝑠𝑡𝑚 ), (1)
𝑑𝑛 𝑑𝑡 = 1 𝜏 (𝑔(𝑣) -𝑛), (2)
𝑓(𝑣) = { 𝑎 𝑓𝑛 (𝑣 -𝑏 𝑓𝑛 ) 2 + 𝑐 𝑓𝑛 (𝑣 < 0) 𝑎 𝑓𝑝 (𝑣 -𝑏 𝑓𝑝 ) 2 + 𝑐 𝑓𝑝 (𝑣 ≥ 0) , (3)
𝑔(𝑣) = { 𝑎 𝑔𝑛 (𝑣 -𝑏 𝑔𝑛 ) 2 + 𝑐 𝑔𝑛 (𝑣 < 𝑟 𝑔 )
𝑎 𝑔𝑛 (𝑣 -𝑏 𝑔𝑛 )

2 + 𝑐 𝑔𝑛 (𝑣 ≥ 𝑟 𝑔 ) .

(

) 4
Here 𝑣 represents the membrane potential, 𝑛 is a variable that reflects the activities of hyperpolarizing ionic channels. 𝑎 𝑥𝑦 , 𝑏 𝑥𝑦 , and 𝑐 𝑥𝑦 ( 𝑥 = 𝑓, 𝑔 and 𝑦 = 𝑛, 𝑝) are parameters. The Parameter 𝐼 0 is a bias constant and 𝐼 𝑠𝑡𝑚 represents the input stimulus.

The only nonlinearity in the DSSN model is a quadratic function. Thus, solving this model using Euler's method needs one multiplication operation per step if the parameters are carefully selected [START_REF] Li | An FPGA-Based Silicon Neuronal Network with Selectable Excitability Silicon Neurons[END_REF][START_REF] Nanami | Simple Cortical and Thalamic Neuron Models for Digital Arithmetic Circuit Implementation[END_REF] . As multiplication operation requires relatively large resources in digital circuits, the DSSN model is suitable for digital silicon neuronal networks.

Synaptic model

Postsynaptic current (PSC) is a current inserted to the postsynaptic cell that induces temporal change in its membrane potential. The PSCs generated by a pulse stimulus to chemical synapses are able to be modeled by the alpha function with double-exponential generalization [START_REF] Amari | The Handbook of Brain Theory and Neural Networks[END_REF] . In our network model, the PSP model was simplified as follows.

𝑑𝑥 𝑑𝑡 = -𝑥 𝜏 .

(

) 5
Here 𝑥 represents PSC generated by a chemical synaptic transmission. It is reset to 𝑤, the connection strength of the synapse, when the membrane potential of the presynaptic neuron exceeds 0. The initial value of 𝑤 is 6x10 -10 . The PSCs in the synaptic connections to a postsynaptic are summed up as follows.

𝐼 𝑠𝑡𝑚 = ∑ 𝑥 𝑖 . (6)
Here i represents index of presynaptic neurons.

STDP algorithm and architecture

Synaptic plasticity means the magnitude of synaptic efficacy to change in response to the activities of pre and postsynaptic neurons.

The STDP is a biological rule that adjusts the strength of synaptic connections (𝑤) based on the relative timing of the spikes in a pre and postsynaptic neurons. Recently it has been shown how STDP rule play a key role by detecting repeating patterns and generating selective response to them [START_REF] Grassia | Spike pattern recognition using artificial neuron and spike-timingdependent plasticity implemented on a multi-core embedded platform[END_REF] .

The STDP rule is a most common form of learning rules used in SNNs. Here is the standard exponential STDP equations.

𝑤 𝑛𝑒𝑤 = 𝑤 + ∆𝑤,

∆𝑤 = { 𝑎 + • exp ( 𝑡 𝑗 -𝑡 𝑖 𝜏 + ) , 𝑖𝑓 𝑡 𝑗 ≤ 𝑡 𝑖 (𝐿𝑇𝑃) -𝑎 -• exp (- 𝑡 𝑗 -𝑡 𝑖 𝜏 -) , 𝑖𝑓 𝑡 𝑗 > 𝑡 𝑖 (𝐿𝑇𝐷) . ( (7)
) 8
Here ∆𝑤 is the modification of the synaptic weight.

When a post-synaptic spike arises after a pre-synaptic spike (𝑡 𝑗 ≤ 𝑡 𝑖 ), the connection is reinforced (long-term potentiation (LTP), ∆𝑤 > 0), whereas in the opposite case it is weakened (long-term depression (LTD)).

Implementation

The overall architecture of our SNN is shown in Fig. 1.The DSSN, STDP, and PSC blocks were implemented on a FGPA chip. In this chip, single postsynaptic neuron is connected with 500 input afferent neurons. Stimulus spike trains were generated by Brian2 on PC and sent to FPGA via Ethernet connection which are explained in 3.4 and 3.5.

Implementation of DSSN

The DSSN model's differential equations are solved by Euler's method and the value of dt is 0.1s. The solver circuit has time division multiplexing (TDM) and pipelined architecture.

Implementation of PSC

As introduced in 2.2, we use a simplified exponential decay to approximate the PSCs. The connection between a postsynaptic neuron and the 500 pre-synaptic neurons was calculated efficiently by using TDM and two-stage pipeline as in the DSSN solver circuit. Operations for the 500 sets of PSCs were finished in 1000 clock cycles, by using single multiplier and single adder. The time step for the summation of PSCs is within 0.1ms if the system clock is faster than 100MHz.

Implementation of STDP learning

Based on the STDP rule, each presynaptic (postsynaptic) spike, induces an LTD (LTP) where the synaptic weight is updated according to Eq. ( 7). An example waveform of PSC is shown in Fig. 3, where the solid and dotted red curve represent the presynaptic and postsynaptic spikes, as well as the black curve is the strength of synaptic connection. For simulating 500 afferents, 500 clock cycles were consumed. Calculation of exponential function and update of synaptic efficacy are executed by TDM.

Ethernet on FPGA

Ethernet is a computer networking technology commonly used in local area networks (LAN). Devices equipped with Ethernet interface are connectable to LAN as well as the internet by supporting common communication protocol, for instance, TCP/IP and UDP. For neuronal network implementation on FPGAs, transmission of the spikes is always an issue particularly when the number of neurons is large. By implementing Ethernet interface and a full hardware protocol stack including IP, UDP and ARP protocol, our SNN on FPGA is capable of receiving input spike patterns from a PC by Ethernet connection and sending report frame back to the PC for monitoring the SNN working status at any time. The configuration and parameter setting are also possible by the Ethernet connection, which contribute to take advantage of FPGA in flexibility and portability.

Pattern recognition

Spatiotemporal spike pattern recognition tasks proposed by Masquelier [START_REF] Masquelier | Competitive STDP-based spike pattern learning[END_REF] were performed on our network. The stimulus pattern required in them were generated by superimposing a spatiotemporal spike pattern at many time points on a background random spike pattern. The former spike patterns (ones to be detected) were generated by PoissonGroup Function in Brian2 (a Python library) on PCs. The background random spike patterns were generated using the same library. The stimulus patterns were stored in files and sent from the PC to the FPGA via the Ethernet connection. An example of learning process is shown in Fig. 6. The postsynaptic neuron started to detect the superimposed patterns at 11s.Our circuit performed this task in real-time. FPGA's system clock was 100MHz.

Conclusion

In this paper, a SNN on an FPGA with STDP learning capability was reported. Stimulus spikes were conveyed in real-time to the FPGA chip via Ethernet connection. It was proved that our SNN is capable of spatiotemporal spike pattern recognition in real-time. In this work, the input patterns were generated on a PC. Perspectives are to generate them in real-time with another SNN. It is also planned to develop an FPGA-FPGA connection bus to expand the scalability of the FPGA-based SNNs, as the scale of the network on an FPGA chip is limited by the amount of on-chip memory used for storing the strength of synaptic connections. An SNN platform with high scalability is advantageous compared to software simulation from the viewpoint of real-time operation and power consumption.

Fig. 1 . 2 .

 12 Fig. 1. Overall Architecture of STDP learning. Fig. 2. Stimulation of a DSSN neuron by a constant input current (0.308).

Fig. 3 .

 3 Fig.3. Waveform of STDP method implementation.

© The 2020 International Conference on Artificial Life and Robotics (ICAROB2020), Feb. 20-23, B-Con Plaza, Beppu, Oita, Japan

Acknowledgements

This work was supported by JSPS KAKENHI Grant Number 15KK0003.


Jason K Eshraghian

Seungbum Baek

Timothée Levi

Takashi Kohno

Said Al-Sarawi

Philip H W Leong

Kyoungrok Cho

Derek Abbott

Omid Kavehei

Nonlinear retinal response modeling for future neuromorphic instrumentation

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Nonlinear Retinal Response Modeling for Future Neuromorphic Instrumentation

Jason K. Eshraghian, Seungbum Baek, Timothée Levi, Takashi Kohno, Said Al-Sarawi, Philip

Leong, Kyoungrok Cho, Derek Abbott and Omid Kavehei

The development of a bio-inspired image sensor that can match the functionality of the vertebrate retina in terms of image resolution, power dissipation and dynamic resolution demands an extremely challenging standard in terms of efficiency and performance [START_REF] Mcintosh | Deep learning models of the retinal response to natural scenes[END_REF][START_REF] Saha | Pervasive computing: A paradigm for the 21st century[END_REF][START_REF] Mead | Neuromorphic electronic systems[END_REF].

At present, functional biological image processing possesses extremely high temporal resolution by asynchronously optimizing the sampling rate of a scene but suffer from low spatial resolution, noise, and an inability to process low frequency content [START_REF] Lichtsteiner | A 128×128 120 dB 15μs Latency Asynchronous Temporal Contrast Vision Sensor[END_REF]. Regardless, the use of dynamic vision sensors that emulate biological retina behaviors are important in performing classification tasks using high-speed motion detection at low-power requirements [START_REF] Linares-Barranco | Dynamic Vision Sensor integration on FPGA-based CNN accelerators for high-speed visual classification[END_REF][START_REF] Moini | An Insect Vision-Based Motion Detection Chip[END_REF][START_REF] Naeini | A Novel Dynamic-Vision-Based Approach for Tactile Sensing Applications[END_REF]. The difficulty in pushing image sensor performance to meet the specifications of the retina lies in the biological complexities of the underlying neural networks. Therefore, a prerequisite for hardware mapping of biological vision systems is understanding the bioprocessing between various retina cells. Here, we present a biologically plausible cellular network retina simulator where light stimuli are accepted at the front-end array of photoreceptor cells, and pass through a system of nonlinear integral equations derived from experimental voltage-and current-clamp data used to calculate the response of each sequential retina cell. It was shown in [START_REF] Eshraghian | Formulation and implementation of nonlinear integral equations to model neural dynamics within the vertebrate retina[END_REF] that integrating the system improves performance when compared with conventional numerical solvers. We expect that by placing our simulator in the hands of interdisciplinary researchers, especially computational neuroscientists, mathematical scientists and machine learning practitioners, it will foster the development of more efficient representations of visual inputs in deep learning applications.

Retinomorphic Processing

The retina can be approximated into five levels of cellular layers through which incident light passes: the rod and cone photoreceptors (outer boundary of the retina), horizontal cells, bipolar cells, amacrine cells and ganglion cells (inner boundary). When light enters the eye, it encounters two types of photoreceptors: the rod cells (highly sensitive in photopic conditions and that enable color vision) which convert light into nerve signals [START_REF] Kefalov | Role of visual pigment properties in rod and cone phototransduction[END_REF] and the cone cells (sensitive in scotopic conitions and assist with vision in dimly lit scenes). In retina modeling the electromagnetic frequency of a photon is often ignored due to the significant complexity associated with measurement of spectral sensitivity in a photoreceptor mosaic, in conjunction with simulating a large-scale cellular network. In the RGB color space, this would triple the throughput of a system. However, with the development of techniques that significantly reduce the computational resources needed for spike-train generation such as pre-integration of compartmental retina models, we can demonstrate a real-time retinal mode that reduces simulation runtimes by up to 50% over equivalent state-of-the-art methods without compromising signal accuracy in using larger time steps. Notably, this model not only produces ganglionic cell spike-trains that are dependent upon the intensity and wavelength of incoming light, but also the intermediary action potentials and graded potentials produced by the electrical and chemical intra-and inter-cellular interactions.

Biologically plausible models of image capture and processing seek to bridge the physical gap between sensing and computation. In doing so, advantages in comparison to traditional cameraprocessor approaches manifest in the form of heightened parallelism, small chip size and low power dissipation, and high dynamic range [START_REF] Eshraghian | Neuromorphic Vision Hybrid RRAM-CMOS Architecture[END_REF], [START_REF] Vanarse | Neuromorphic Engineering -A Paradigm Shift for Future IM Technologies[END_REF]. To achieve this standard of neuromorphic hardware implementation, there are several barriers to overcome. These include the lack of access to experimental data from the vertebrate retina, the difficulty of processing large-scale compartmental models without introducing simplified assumptions in real-time, and the limited interconnectivity in silicon across a 2D planar chip.

In this work, we seek to address the issue of real-time processing without compromising accuracy. We provide open-access to a retinal simulation platform for researchers at the end of this paper, to facilitate a wholly integrated real-time model of a measurable retinal network with adaptive parameters, including a variation of cell membrane capacitance, maximum conductance, reversal potential, ion channel selection, and electrical/chemical synaptic parameterization. We seek to open up the frontiers for accessible and reproducible benchmarks in neuromorphic instrumentation for vision processing, in both individual and populations of neuronal measurements to provide a more in-depth insight into the role of the various physical, chemical and electrical properties of the retinal cellular network. We will first demonstrate the high-level architecture of our retinal measurement platform as depicted in Fig. 1 before looking at the constituent low-level compartmental blocks, and the resultant data procured from the retinal network.

Signal Processing and Computation Architecture

With higher precision and greater parameterization comes slower processing. Computer vision algorithms are often implemented without regard to the computer architecture or hardware they are designed to run on, and so similar structures and principles should ideally exist between algorithm and hardware. With the construction of a large system of integral equations it becomes necessary to accelerate computational efficiency, where from a mathematical perspective, a system of ordinary differential equations (ODEs) and integral equations describe identical phenomena. However, from an algorithmic-hardware co-design analysis, three salient benefits arise:  Ability to solve the large system of nonlinear integral equations in real-time using small time-steps of sub-nanosecond order.

 Reduced reliance on off-chip memory, as hardware implementation of derivatives require memory and a delay line, whereas our integral formulation only requires summations for correct implementation.

 Biophysically accurate representation of the retina at a cellular level.

The formulation of a conductance-based compartmental retinal model initiates at the first input layer of photoreceptor cells, passes through three hidden layers, and triggers the output ganglion cells which respond with a spike-train sent to the visual cortex for further processing and interpretation. The system of ODEs can be generically characterized by:

𝑑[Ca 2+ ] 𝑑𝑡 = 𝐼 Ca (𝑉) + 𝐺(𝑉; [Ca 2+ ]), (1)
𝑑𝜒 𝑑𝑡 = (1 -𝜒)𝛼 𝜒 (𝑉) -𝛽 𝜒 (𝑉)𝜒, (2)
where 𝐶 𝑚 is the membrane capacitance, V is the membrane potential, 𝐼 syn (𝑉) = 𝑔(t)(𝐸 syn -𝑉) is the synaptic current input, 𝐸 syn is the reversal potential, and 𝐼 ion is the sum of all ionic voltage and calcium-dependent ionic currents. Here, [Ca 2+ ] is the intracellular concentration of free calcium, 𝐼 𝐶𝑎 (𝑉) denotes all voltage-activated calcium currents, while G represents all other voltage and calcium-dependent cellular processes. The gating variable is denoted χ, and forward-backwards propagation rates are represented by 𝛼 𝜒 and 𝛽 𝜒 , respectively. While the complete and detailed system of equations, initial conditions, biasing and physical descriptions are subsumed within the high-level description, their details can be found in an online appendix with a link provided in the references [START_REF] Baek | Mathematical Modeling of Retinal Cells[END_REF].

The signal processing architecture of the integrated form of equations ( 1) -( 3) is shown in Fig. 1. The architecture provides a visual understanding of how the system of nonlinear integral equations behaves within a nonlinear time-invariant system (that is, the light input is the timevarying component measured in the number of rhodopsin, Rh* in Fig. 2, and photoisomerization molecules activated per second).

The computational mapping of a photocurrent signal from the input photoreceptor layer through to the ganglion cell output layer requires an input conductance g(t) that is weighted by the spectral sensitivity of the cone cell it comes into contact with. The resulting signal is passed to the primary voltage/ionic current evolution block where it is used to derive the synaptic current.

The summation of currents and gating variables is fed to the total calcium current in a secondary block that describes the evolution of calcium concentration. The signals are combined and integrated to generate new values for membrane potential and calcium at the next time step ti+1.

This formulates the high-level architecture for the retinal network. A sequential demonstration of the first two layers of photoreceptor (rod and cone) and bipolar cells is provided in Fig. 2.

Retinal Model

Spectral Sensitivity

Light incident on the retina is a mixture of different wavelengths, and at least three different subclasses of cone photoreceptor cells have been experimentally verified to be wavelengthdependent. Each of these subclasses can be classified based on the wavelength sensitivity of the photopigment in its outer segment: long-wavelength (L), middle-wavelength (M) and shortwavelength (S) peak sensitive cone cells.

The photoreceptor mosaic of varying spectral sensitivities is implemented in the proposed model as input dependent weighting coefficients, where normalized sensitivity is represented by weights ŜS, ŜM and ŜL that are approximate models of cone cell responses across the spectrum of visible wavelength, based on corneal measurements presented in [START_REF] Stockman | Spectral sensitivitesof the human cones[END_REF]. These models were fit using least squares regression with the 5 th order polynomial, which quantitatively match the precision of the rest of the model:

Ŝ 𝑆 = { -1.
3 × 10 -9 + 3.1 × 10 -6 𝜆 4 -2.9 × 10 -3 𝜆 3 +1.3𝜆 2 -310𝜆 + 2.9 × 10 4 , 400 ≤ 𝜆 ≤ 530 0, otherwise Ŝ 𝑀 = { -0.099𝑧 5 + 0.24𝑧 5 -0.4𝑧 3 + 0.94𝑧 2 -0.37𝑧 + 0.87, 400 ≤ 𝜆 ≤ 530 0, otherwise

Ŝ 𝐿 = { 0.12𝑦 5 + 0.21𝑦 4 -0.55𝑦 3 -0.88𝑦 2 +0.59𝑦 + 0.9, 450 ≤ 𝜆 ≤ 620 0, otherwise

where 𝜆 is wavelength in nanometres, 𝑧 = 𝜆-5.4×10 2

76

and 𝑦 = 𝜆-5.5×10 2

85

. Another reason for using polynomial forms, more generally, is the simplicity of using artificial neural networks to encode polynomial functions, which becomes essential for deep learning applications and predictive retinal networks.

The cone cell variations in spectral sensitivity demonstrate that the biological retina does not respond equally to all types of light, and does not have a uniform sensitivity across the visible spectrum. In addition to this, there is further variation across the photoreceptor mosaic (our field of view) to light intensity where the fovea is considered to be the centre. Therefore, the Westheimer linespread function can be used to capture the characteristics of polychromatic light more accurately:

Ŝ = 0.47𝑒 -3.3𝑖 2 + 0.53𝑒 -0.93|𝑖| , (7)
which now introduces a dependence on not only photon wavelength but also on the position upon the photoreceptor array where the light impinges, where i is the position on the retina (on a 2-dimensional axis) specified in terms of minutes of visual angle from the fovea. We can calculate the different weights based on input wavelength 𝜆 and position on which light falls on the retina i for a short-wavelength peak sensitive cone cell by taking the root of ( 4), and multiplying it by [START_REF] Naeini | A Novel Dynamic-Vision-Based Approach for Tactile Sensing Applications[END_REF], which has been graphically illustrated in Fig. 3. This is a near perfect fit of the benchmark data of spectral sensitivities across the photoreceptor mosaic in [START_REF] Wandell | Foundations of Vision[END_REF], and also demonstrates that peak sensitivity occurs at i = 0 for a wavelength of 𝜆 = 440 nm. To qualitatively summarize the function of ( 4) -( 7), the input light is weighted using our spectral sensitivity model based on the following three features:

 whether light targets a L, M or S cone cell;

 the wavelength of the incoming light, and  the visual angle of incoming light concerning the targeted cell.

A Conductance-Based Compartmental Model of the Retina

Thus far, we have elucidated the time-invariant architecture of our model, the wavelength and spatial-dependence on impinging light. The computational bottleneck of most retinal systems, however, lies in emulation of the high-dimensional spaces that much of the molecular activity within biological cells occurs at. With neural modeling, there are often occasions where, due to an absence of data for a specific cell in a species, one is required to adopt and build neural models whose ionic mechanisms are not from the same class, but from related organisms. The reality is that no complete dataset that wholly characterizes the retinal cellular network of a single animal exists. In light of this, we have procured a substantial collation of experimental data that describes ionic channels, membrane channel densities, and temporal synaptic dynamics for each cell in the retina and performed our own fitting on information that was not previously characterized, in our endeavour to build a unified model with correct biophysical properties for the vertebrate retina. Table 1 depicts the sources of all experimental data that make up the system of nonlinear ODEs from ( 1) -( 3), which is subsequently computed as a system of a total of 272 nonlinear integral equations and parameters in Fig. 1. Some of these sources have utilized their own collations of various species in developing either their models or experimental data, in which case, these references have been grouped in accordance with color. Notably, the squid is the only species on the list which is not classified as posessing a vertebrate retina. However, the Hodgkin-Huxley conductance-based model, while derived from the squid, can be extended to the vertebrate retina due to the likeness in the ionic mechanisms underlying the initiation and propagation of action potentials [START_REF] Hodgkin | A quantitative description of membrane current and its application to conduction and excitation in nerve[END_REF]. The parameterizable features of electrical and chemical synapses between neurons include the maximum conductance, the reversal potential, the slope in millivolts, and the time constant. The detailed references from Table 1 can be found in the appendix in addition to all equations and variables [START_REF] Baek | Mathematical Modeling of Retinal Cells[END_REF].

Simulation Platform

The simulation platform was developed in the C language, the user interface in Microsoft Foundation Class, and graph plotting using Kst and is recommended for use on Windows 10.

Full-scale network simulation is depicted in Fig. 4 with decoupled rod and cone cell pathways, and the coupled pathway which exhibits horizontal lateral connections through electrical synapses, in response to a spatially uniform pulse of light that is activated for a duration of 20 ms. Importantly, the simulation was conducted in real-time at a biological time scale, and all listed measurements in Fig. 4 present a new method to attain benchmarks for vertebrate retina signal extraction under a flexible user-configurable input. The light input is first passed through a weighted conditional switch in Fig. 1 as determined by the wavelength of light input, and the type of cone cell. The polynomial/linespread response (demonstrated in Fig. 3 under a specific test case) is processed using the integral form of (1) -

, through repeated cycles of iterating time steps until the solution converges. Table 2 provides a comparative analysis of speed between various numerical solvers, namely the classical fourth-order Runge-Kutta (RK4) method, the more atypical ODE15s solver, and the midpoint methods the latter two being our points of comparison for a conventional system of ODEs. In addition to collectively providing the fundamental basis for our retinal model, the studies undertaken in Table 2 are used as a basis of comparison for a qualitative verification of accuracy, while a quantitative analysis is undertaken in Fig. 5. To quantify the accuracy of the proposed model, the benchmark used was the conventional ODE15s solver from the MATLAB ODE suite at a very small time step of 1.00×10 -9 s. The error for the ganglion membrane potential and the calcium concentration accumulates over time, which arises due to slight misalignments between spike patterns leading to an overstatement of absolute error. When the data is manually timescaled such that peaks of spikes match, the error is bounded and ceases to accumulate.

A range of other conventional ODE solvers were tested (namely, ODE45, ODE23s and ODE113). ODE45 and ODE113 are well-suited for solving non-stiff systems and require unreasonably low time steps for stable and accurate solutions. This significantly increased computation time. ODE15s and ODE 23s are both appropriate for stiff systems, but ODE23s tends to only operate effectively with crude error tolerances. As such, ODE15s was selected as both fastest and most accurate of all the solvers, and we selected only the best performing method as the benchmark for comparison against our proposed system of nonlinear integral equations. Ultimately, the RK4 method is consistently faster than the ODE solver for fixed time steps demonstrating the feasibility of generating real-time predictive retinal measurements.

The values in Table 2 have been normalized to the value of the ODE15s execution time with a time step of 1.00×10 -1 s. One discrepancy in the results has been highlighted and occurs with the midpoint method which is computationally slower at a time step of 1.00×10 -1 s than 1.00×10 - Fig. 5 Absolute error over biological time for the S-cone cell pathway. In all cases, cells display a maximum error at the same time the stimulus is applied for 20 ms.

2 s. This is explained as a result of the time step simply being too large at 1.00×10 -1 s to characterize the full spike at the spike-train output of the ganglion cell, which occurs on the order of milliseconds, and can be verified by reference to the ganglion cell simulation shown by curve (XV) in Fig. 4, which depicts both decoupled and coupled retinal pathways of the ganglion membrane potential in millivolts.

Conclusions

A biologically plausible integral based simulation platform is described in this paper. The simulator is expected to both provide a benchmark for computational neuroscientists who need fast and accurate predictive networks of the retina, and to also assist fundamental studies in mathematics and biomedical sciences that are exploring how the retina is processing information to inspire new ways to introduce sparsity into computer vision. While we are still some distance from a perfect one-to-one model of the retina, we anticipate this simulation platform will help accelerate those studies in order to gain a more in-depth insight into the roles of constituent cells and their associated parameters, and how they may leverage the use of spiking neural networks in assisting other maturing deep learning tasks. Longer-term research motivations can make use of biological real-time simulation for use in bio-hybrid experimentation. The simulator is available for download at the following link: http://ccns.cbnu.ac.kr/paper/1611/retinamodelsimulator.zip.

Fig. 1

 1 Fig. 1 Higher-level architecture of photocurrent signal processing and computation.

Fig. 2

 2 Fig.2Front-end simulation results of rod cell, RGB (LMS) cone cells, and bipolar for a given light stimulus. Here, g(t) is the input conductance of the bipolar cell ion channel.

Fig. 3

 3 Fig. 3 Simulation of normalized spectral sensitivity as a function of light wavelength and the visual angle away from the targeted S-cone cell as per the proposed model.

Fig. 4

 4 Fig.[START_REF] Lichtsteiner | A 128×128 120 dB 15μs Latency Asynchronous Temporal Contrast Vision Sensor[END_REF] Photoreceptor stimulation to a spatially uniform pulse of light for a duration of 20 ms, shown across varying timescales. This is the equivalent of 100 Rhodopsin molecules activated per second (Rh*/sec) for rod cells, and 50,000 photoisomerizations per cone cell per second. For simplicity, only the S-cone cell was triggered. The ganglion membrane potential graph in (XV) includes an overlay of spike trains produced by intact cell recordings (current clamp) in the tiger salamander retina[START_REF] Fohlmeister | Modeling the repetitive firing of retinal ganglion cells[END_REF]. Note that amacrine cell simulations are shown in the rod cell pathway in plots VI and VII. Lateral connections via electrical synapses are modeled between the input photoreceptors. Horizontal cells facilitate signal transmission across different pathways, though have been omitted from these simulations to focus on the behavior of a pathway of a single event.

Table 1 -Feature collation for retinal modeling

 1

	Species	Function	Literature
	Collated	Bandpass filtering of rod photoreceptor network	Kamiyama et al. (2009)
	Collated	Ionic current in rod photoreceptor network	Kamiyama et al. (1996)
	Salamander	Voltage-and calcium-activated current of rod cell inner	Bader et al. (1982)
		segment
	Lizard	Calcium and calcium-dependent chloride currents in cone	Maricq et al. (1988)
		photoreceptor cells
	Salamander	Ionic channels of cone photoreceptors	Barnes et al. (1989)
	Toad	Phototransduction in rod cells	Torre et al. (1990)
	Salamander	The conductance of rod cells	Baylor et al. (1986)
	Newt	Phototransduction in rod cells	Forti et al. (1989)
	Turtle	Electrical responses of cone cells	Baylor et al. (1974)
	Collated	Ionic current in bipolar cells	Usui et al. (1996)
	Goldfish	Hyperpolarization of activated current in the cell body	Kaneko et al. (1985)
	Goldfish	Calcium current in the axon terminal	Tachibana et al. (1991)
	Goldfish	Ca-dependent chloride current in the axon terminal	Tachibana et al. (1993)
	White bass	Delayed rectifying potassium current in the cell body	Lasater (1988)
	Axolotl	GABA-induced current in the axon terminal	Attwell et al. (1987)
	Axolotl	Glutamate-induced current in the dendrite	Attwell et al. (1987)
	Salamander	Glutamate-induced current in the dendrite	Nawy et al. (1990)
	Dogfish	Glutamate-induced current in the dendrite	Shiells et al.

Table 2 -Execution Time as a Function of Time Steps Time step (s)

 2

			Normalized Execution Time
			S-Cone Pathway
		Midpoint	RK4	ODE15s
	1.00×10 -1	2.33	1.00	1.00
	1.00×10 -2	1.67	2.67	168
	1.00×10 -3	16.0	20.0	183.3
	1.00×10 -4	180.7	231.0	375.3
	1.00×10 -5	1,182	1,738	2,139
	1.00×10 -6	11,150	17,194	21,184
			Normalized Execution Time
	Time step (s)		Rod Pathway
		Midpoint	RK4	ODE15s
	1.00×10 -1	3.25	1.00	1.00
	1.00×10 -2	1.5	2.00	1,100
	1.00×10 -3	13.3	16.5	1,190
	1.00×10 -4	139.8	186.3	1,268
	1.00×10 -5	1,044	1,518	2,038
	1.00×10 -6	10,049	14,899	15,108

Acknowledgments

This work was supported by the National Research Foundation of Korea Grant funded by the Korean Government under Grant NRG-2017R1D1A1A09000613, by the DFAT Australia-Korea Foundation under Grant AKF00640, and by iDataMap Corporation Pty. Ltd. Australia.


Claudio Attaccalite

Maurizia Palummo

Elena Cannuccia

Myrta Grüning

How strong is the Second Harmonic Generation in single-layer monochalcogenides? A response from first-principles real-time simulations

published or not. The documents may come

How strong is the second harmonic generation in single-layer monochalcogenides: A response from first-principles real-time simulations

Claudio Attaccalite, Maurizia Palummo, Elena Cannuccia, Myrta

Second Harmonic Generation (SHG) of single-layer monochalcogenides, such as GaSe and InSe, has been recently reported [2D Mater. 5 (2018) 025019; J. Am. Chem. Soc. 2015, 137, 79947997] to be extremely strong with respect to bulk and multilayer forms. To clarify the origin of this strong SHG signal, we perform first-principles real-time simulations of linear and non-linear optical properties of these two-dimensional semiconducting materials. The simulations, based on ab-initio many-body theory, accurately treat the electron-hole correlation and capture excitonic effects that are deemed important to correctly predict the optical properties of such systems. We find indeed that, as observed for other 2D systems, the SHG intensity is redistributed at excitonic resonances. The obtained theoretical SHG intensity is an order of magnitude smaller than that reported at the experimental level. This result is in substantial agreement with previously published simulations which neglected the electron-hole correlation, demonstrating that many-body interactions are not at the origin of the strong SHG measured. We then show that the experimental data can be reconciled with the theoretical prediction when a single layer model, rather than a bulk one, is used to extract the SHG coefficient from the experimental data.

I. INTRODUCTION

Nowadays there is a considerable interest in the excited state properties of 2D materials which can be distinctly different from that of their bulk counterpart. [START_REF] Autere | Nonlinear optics with 2d layered materials[END_REF] Strongly bound excitons often characterize the optical spectra of two-dimensional (2D) semiconducting crystals and in several cases, intense nonlinear optical spectra have been observed [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF][START_REF] Säynätjoki | Ultra-strong nonlinear optical processes and trigonal warping in mos 2 layers[END_REF] such as strong second-harmonic generation [START_REF] Malard | Observation of intense second harmonic generation from mos2 atomic crystals[END_REF][START_REF] Kumar | Second harmonic microscopy of monolayer mos2[END_REF][START_REF] Janisch | Extraordinary second harmonic generation in tungsten disulfide monolayers[END_REF][START_REF] Ribeiro-Soares | Second harmonic generation in wse 2[END_REF][START_REF] Zhou | Strong second-harmonic generation in atomic layered gase[END_REF] (SHG)-typically up to one order of magnitude larger than in conventional nonlinear crystals. For this reason these materials are of potential technological interest as frequency converters in nanophotonic circuits. [START_REF] Majumdar | Hybrid 2d material nanophotonics: A scalable platform for low-power nonlinear and quantum optics[END_REF] Furthermore, the SHG is an excellent spectroscopic tool for the imaging and characterization of 2D materials [START_REF] Säynätjoki | Ultra-strong nonlinear optical processes and trigonal warping in mos 2 layers[END_REF][START_REF] Malard | Observation of intense second harmonic generation from mos2 atomic crystals[END_REF][START_REF] Kumar | Second harmonic microscopy of monolayer mos2[END_REF][START_REF] Zhou | Strong second-harmonic generation in atomic layered gase[END_REF][START_REF] Li | Probing symmetry properties of few-layer mos2 and h-bn by optical second-harmonic generation[END_REF][START_REF] Hsu | Second harmonic generation from artificially stacked transition metal dichalcogenide twisted bilayers[END_REF][START_REF] Yin | Edge nonlinear optics on a MoS atomic monolayer[END_REF] and understand the origin of the signal from theoretical point of view is of fundamental importance for the interpretation.

Recently, extremely large values for the SHG of singlelayer monochalcogenides have been reported: for InSe 13 a value of 6.39 nm/V at 1.55 eV and for GaSe 8 a value of 2.4 nm/V at 1.02 eV. To put these values into perspective, the SHG of single-layer MoS 2 is of the order of 0.1-0.4 nm/V 4,14 at resonance, that of bi-and tri-and multilayer [START_REF] Jie | Layer-dependent nonlinear optical properties and stability of non-centrosymmetric modification in few-layer gase sheets[END_REF][START_REF] Tang | Layer-and frequency-dependent second harmonic generation in reflection from gase atomic crystals[END_REF][START_REF] Karvonen | Investigation of second-and thirdharmonic generation in few-layer gallium selenide by multiphoton microscopy[END_REF] GaSe ranges from 9-90 pm/V, and that of bulk GaSe is about 60 pm/V. A similar large SHG of 10 nm/V was reported for WS 2 and WSe 2 monolayers. [START_REF] Janisch | Extraordinary second harmonic generation in tungsten disulfide monolayers[END_REF][START_REF] Ribeiro-Soares | Second harmonic generation in wse 2[END_REF] However, the extremely large SHG coefficients for single-layer monochalcogenides are not confirmed by first-principle calculations [START_REF] Lin | Indium selenide monolayer: a two-dimensional material with strong second harmonic generation[END_REF][START_REF] Hu | Layerindependent and layer-dependent nonlinear optical properties of two-dimensional gax (x= s, se, te) nanosheets[END_REF] which give consistently values smaller by one order of magnitude.

The discrepancy between the theoretical and exper-imental results may be due to many-body effects that were neglected in the computational studies of Refs. 18 and 19. In fact, electron-electron interaction and excitonic effects are known to play a crucial role in the description of optical properties, especially in the case of low-dimensional materials. Indeed the incomplete screening of the electron-electron interaction and the quantum confinement may further enhance these effects with respect to bulk. Specifically for single-layer monochalcogenides, the formation of strong bound excitons have been recently shown in Ref. 20. Further, their peculiar optical properties have been attributed to saddle-points excitations [START_REF] Luo | Saddle-point excitons and their extraordinary light absorption in 2d β-phase groupiv monochalcogenides[END_REF] that originate from the Mexican hat-like band structure near the Γ point. This work addresses the question of whether the extremely large SHG of single-layer monochalcogenides is due to strongly bound excitons. To this end, by using a first-principle real-time approach based on many-body theory (Sec. II), we calculate the electronic structure, di-arXiv:1903.06563v2 [cond-mat.mtrl-sci] 24 Jul 2019 electric function (Sec. III A) and the SHG (Sec. III B) of InSe and GaSe monolayers and, in order to single out many-body effects, we compare with the results at the independent-particle level. Previous studies on h-BN, MoS 2 and other semiconducting 2D materials found a significant enhancement at resonance of the SHG signal due to excitonic effects-up to 2-4 times depending on the system. [START_REF] Grüning | Second harmonic generation in h-bn and mos2 monolayers: Role of electron-hole interaction[END_REF][START_REF] Attaccalite | Strong second harmonic generation in sic, zno, gan twodimensional hexagonal crystals from first-principles manybody calculations[END_REF] For single-layer monochalcogenides, we find that intensity is only slightly redistributed by many-body effects respect to the independent-particle approximation. We conclude (Sec. IV) that many-body effects do not account for the difference between theoretical predictions [START_REF] Lin | Indium selenide monolayer: a two-dimensional material with strong second harmonic generation[END_REF][START_REF] Hu | Layerindependent and layer-dependent nonlinear optical properties of two-dimensional gax (x= s, se, te) nanosheets[END_REF] and experimental estimates. [START_REF] Zhou | Strong second-harmonic generation in atomic layered gase[END_REF][START_REF] Zhou | Inse monolayer: synthesis, structure and ultrahigh second-harmonic generation[END_REF] We turn then the attention to the experimental results and argue that the extremely large SHG coefficient is an artifact of the model assumed to extract the estimate from the experimental data.

II. THEORETICAL FRAMEWORK AND COMPUTATIONAL DETAILS

The structural optimization and electronic structure have been calculated using density-functional theory (DFT). The crystal structure of the GaSe and InSe monolayers (Fig. 1) has been obtained from the corresponding bulk structure using the same in-plane lattice parameters (a = 0.3743 nm and a = 0.4 nm respectively for the GaSe [START_REF] Schwarz | Effect of pressure on the structural properties and electronic band structure of gase[END_REF] and the InSe monolayer [START_REF] De Blasi | Large inse single crystals grown from stoichiometric and nonstoichiometric melts[END_REF] ). Since experiments on GaSe and InSe were performed on different substrates, the lattice constant of these materials could depend on the substrate interaction. Therefore we decided to keep fixed the lattice constants to the bulk values and optimizing only the atomic positions that are not fixed by symmetries, namely the z coordinates of the atoms. In order to simulate an isolated monolayer, we used a 35 a.u. supercell in the z-direction.

DFT calculations have been performed with the Quantum-Espresso code 28 using the Perdew-Burke-Ernzerhof 29 (PBE) functional and the scalar-relativistic optimized norm-conserving Vanderbilt 30 pseudopotentials from the PseudoDojo repository (v0.4). [START_REF] Van Setten | The pseudodojo: Training and grading a 85 element optimized norm-conserving pseudopotential table[END_REF] The valence configuration for each pseudo-atoms are: 3d 10 4s 2 4p 1 for the Ga, 4d 10 5s 2 5p 1 for the In and 4s 2 4p 4 for Se. We used a shifted 18 × 18 × 1 k-point sampling for the ground-state, a plane-wave cutoff of 90 Ry for the structural optimization and 70 Ry for the band-structure calculations.

Calculations of the quasi-particle energies and optical susceptibilities have been performed using the Yambo code. [START_REF] Sangalli | Many-body perturbation theory calculations using the yambo code[END_REF] The quasiparticle corrections to the fundamental band gap have been calculated from the standard G 0 W 0 approximation 33 with the Godby-Needs plasmonpole model [START_REF] Rw Godby | Metal-insulator transition in kohn-sham theory and quasiparticle theory[END_REF] and applied as a rigid shift to all the bands. We used a 24 × 24 × 1 k-point sampling, a 5 Hartree cutoff for the dielectric function, and 200 total bands for Green's function expansion. We calculate the gap correction at Γ point and then shifted all conduction bands by this amount. In this way the quasiparticle band structure has the same band ordering and band-width of the DFT one, with the only difference that the gap has been corrected within the G 0 W 0 approximation. Convergence with the number of conduction bands was accelerated by means of the Bruneval-Gonze terminator. [START_REF] Bruneval | Accurate gw selfenergies in a plane-wave basis using only a few empty states: Towards large systems[END_REF] Optical absorption spectra have been obtained by the solution of the Bethe-Salpeter equation 36 using a basis of electron-hole pairs for which we considered 6 valence and 7 conduction states for both GaSe and InSe monolayers. We verified that increasing the number of bands in the calculations does not change the spectra in the energy range considered here. Non-linear susceptibilities have been obtained from the real-time evolution of Bloch-electrons in a uniform time-dependent electric field following the approach proposed in Ref. 24 and 37. The effective Hamiltonian for the Bloch-electrons, derived from many-body theory, contains both the electron-hole attraction (through a screened-Coulomb term) and exchange (through a Hartree term) needed to describe excitons and local-field depolarization effects. The whole framework is based on DFT and quasiparticle corrections are included at the level of G 0 W 0 . In the linear-response limit, the approach is equivalent to the solution of the BSE on top of the G 0 W 0 electronic structure [START_REF] Attaccalite | Realtime approach to the optical properties of solids and nanostructures: Time-dependent bethe-salpeter equation[END_REF] and we therefore refer to this level of theory as G 0 W 0 +BSE for both linear and nonlinear susceptibilities. Within this framework, we can exclude specific terms from the effective Hamiltonian to investigate how they affect the SHG. In particular, we con-sider the independent-particle approximation (IPA)-no electron-hole correlation; the random-phase approximation (RPA)-only electron-hole exchange. Further, we perform the simulations on top of either the DFT or the G 0 W 0 electronic structure. For the screened-Coulomb and Hartree terms, we use the same parameters as in the BSE calculations.

After the integration of the equation-of-motion, the nonlinear susceptibilities χ (n) are extracted using Fourier techniques from the macroscopic polarization Pcalculated as a dynamic Berry-phase [START_REF] Souza | Dynamics of berry-phase polarization in time-dependent electric fields[END_REF] -and its expansion in power of the total electric field E:

P i = χ (1) ij E j + χ (2) ijk E j E k + χ (3) ijkl E j E k E l + O(E 4 ) , (1)
where the i, j, k, l subscripts refer to Cartesian components of the field and polarization.

Specifically, to obtain a single SHG spectrum we perform a series of simulations of about 75 fs for a monochromatic electric field with frequency ranging from 0.1 eV to 7.0 eV. We integrate the equation-of-motion using the Crack-Nicholson method [START_REF] Crank | A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type[END_REF] with a time-step of 0.01 fs. Since the turning-on of the electric-field excites all eigenfrequencies of the system, we introduce a phenomenological decoherence non-Hermitian operator corresponding to a Lorentzian broadening of the spectrum of 0.2 eV. 41 The contribution of the eigenfrequencies to the signal decays exponentially in time with the decay constant being the inverse of the dephasing time. After five times the dephasing time (about 33 fs) the eigenfrequency signal is negligible with respect to the SHG signal and we then perform the Fourier analysis as detailed in Ref. 37.

The static limit (ω = 0), corresponding to the outmost left point of Figs. 3,4 has been obtained by extrapolation of the values at small frequencies of the χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (ω).

The external electric field is polarized along the y direction and the polarization is recorded in the same direction, obtaining the χ

(2) yyy = -χ (2) xxy = -χ (2) yyx = -χ (2) xyx = -χ (2)
aab that is the only nonzero component of the secondorder susceptibility tensor for the hexagonal D 3h crystal class, [START_REF] Boyd | Nonlinear Optics[END_REF] being a and b the in-plane crystal axes. To obtain a SHG signal independent of the dimension of the supercell, we rescaled the calculated χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (ω) by an effective thickness of 0.796 nm for the GaSe and 0.832 nm for the InSe-corresponding to half of c lattice parameter of the bulk structures. [START_REF] Boyd | Nonlinear Optics[END_REF]

III. RESULTS

A. Electronic structure and optical properties

In Fig. 2, we report the electronic band structure of both monochalcogenide monolayers and the position of the lowest excitons obtained from the solution of the BSE. Because of the mirror-plane symmetry (z → -z), electronic states near the edges are either even or odd with respect to this symmetry operation. [START_REF] Zólyomi | Electrons and phonons in single layers of hexagonal indium chalcogenides from ab initio calculations[END_REF] As a result, the lowest-energy electron-hole transition-depicted by a white arrow-is optically inactive for in-plane polarized light and only slightly coupled to z-polarized light. [START_REF] Antonius | Orbital symmetry and the optical response of single-layer mx monochalcogenides[END_REF][START_REF] Denis A Bandurin | High electron mobil-ity, quantum hall effect and anomalous optical response in atomically thin inse[END_REF] In the same figure, we also indicate the lowest optical active transition with a black arrow.

The band structure is in good agreement with previous results obtained by Hu et al. [START_REF] Hu | Layerindependent and layer-dependent nonlinear optical properties of two-dimensional gax (x= s, se, te) nanosheets[END_REF] , while we found small differences respect to Antonius et al. [START_REF] Antonius | Orbital symmetry and the optical response of single-layer mx monochalcogenides[END_REF] who reported a larger energy difference between the two top valence bands at the vicinity of Γ. [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (ω = 0), obtained from real-time simulations at different level of approximation: independent-particles approximation (IPA) and the random-phase approximation (RPA) on top of the DFT electronic-structure, the IPA (G0W0) and the real-time BSE (G0W0+BSE) on top of the G0W0 electronic-structure. Units in pm/V.

χ (2) (IPA) χ (2) (RPA) χ (2) (G0W0) χ (2) (
In Fig. 3 and Fig. 4 (bottom panels), we report the imaginary part of the dielectric function of GaSe and InSe for in-plane polarized light. The first excitons, obtained by full diagonalization of the Bethe-Sapeter equations [START_REF] Richard M Martin | Interacting Electrons[END_REF][START_REF] Strinati | Application of the Green's functions method to the study of the optical properties of semiconductors[END_REF] are reported as well. As one can see for both GaSe and InSe we are in presence of bound excitons (i.e. at energies below the fundamental band-gap).The lowest exciton at 3.39 eV in GaSe is dark for light polarized in the plane due to the symmetry of the bands close to the Γ point (while it has a small but non-zero dipole for light polarization in the z direction). The first bright exciton is found at 3.51 eV. Our results are in good agreement with those of Antonius et al. [START_REF] Antonius | Orbital symmetry and the optical response of single-layer mx monochalcogenides[END_REF] For the InSe monolayer (bottom panels of Fig. 4) we obtained similar results. The lowest excitation at 2.47 eV is dark for in-plane polarized light and the first bright exciton is at 2.63 eV. This is in agreement with luminescence measurements that found that the lowest-energy transition is not optically active in InSe monolayer. [START_REF] Denis A Bandurin | High electron mobil-ity, quantum hall effect and anomalous optical response in atomically thin inse[END_REF] Therefore, the lowest excitons do not contribute to the second-order susceptibility for in-plane polarized light.

B. Second Harmonic Generation

In Table I, we report the static χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (0) at different level of theory: IPA, RPA, independent-particle on top of the G 0 W 0 band structure (G 0 W 0 ), and finally the full G 0 W 0 +BSE. At the independent-particle level we found a large static second order susceptibility for both systems. The susceptibility is sensitive to the electronic structure and in particular to the band gap: a larger band gap corresponds to a smaller χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (0) as one can see going from DFT to the G 0 W 0 . Also, the band structure of InSe Arb. Units c) is similar to that of GaSe and the gap at DFT level is 1.82 eV against 2.34 eV of GaSe. As a consequence InSe has a larger χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (0).

² 2 (!) ² 2 (!=2) b) G 0 W 0 + BSE Imag Real Abs 1 2 3 4 5 6 7 Energy [eV] d) ² 2 (!) ² 2 (!=2)
The GaSe result at G 0 W 0 level is compatible with the result of Hu et al. that found a χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (0) = 69.4 pm/V starting from a DFT band structure obtained by means of the Heyd-Scuseria-Ernzerhof hybrid functional that has a gap of 3.34 eV, in-between the PBE 2.34 eV and the G 0 W 0 4.03 eV band gap values. For the InSe monolayer, Lin et al. [START_REF] Lin | Indium selenide monolayer: a two-dimensional material with strong second harmonic generation[END_REF] report a value of 91.3 pm/V, compatible with our result, considering that they assumed a smaller effective thickness. Lin et al. also investigated other possible structures of the InSe monolayer with a stronger SHG but we are not aware of experimental realization of these polymorphs and therefore we limit our analysis to the hexagonal monolayer.

Within the RPA, local field effects-that originate from charge fluctuations induced by the system inhomogeneity-attenuate the second order susceptibility in both GaSe and InSe monolayers. This reduction has been observed as well for other low-dimensional systems and in bulk. [START_REF] Grüning | Second harmonic generation in h-bn and mos2 monolayers: Role of electron-hole interaction[END_REF][START_REF] Attaccalite | Nonlinear optics from an ab initio approach by means of the dynamical berry phase: Application to second-and third-harmonic generation in semiconductors[END_REF] Within the G 0 W 0 +BSE, excitonic effects tend to increase the χ [START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (0) and counteract the reduction from both local-field effects and quasiparticle corrections, in agreement with what observed for bulk semiconductors. [START_REF] Chang | Excitonic effects on optical second-harmonic polarizabilities of semiconductors[END_REF] The SHG intensity at G 0 W 0 +BSE level is for both materials closer to the IPA one.

Figures 3 and4 [(a) and (b) panels] show the imaginary, real part and absolute value of the SHG for GaSe and InSe monolayers within the IPA (a) and G 0 W 0 +BSE (b). They are compared with the imaginary part of the dielectric function at ω and ω/2 (bottom panels) at the same level of theory. From the comparison, the main spectral features in the SHG are easily recognizable as two-photon resonances at the van Hove (IPA) and exciton energies (G 0 W 0 +BSE). We focus here on the features below the onset of the optical absorption, that is in the transparency region (imaginary part of the dielectric function ε 2 ), which is the spectral region of technological interest. For GaSe, comparing the IPA with the G 0 W 0 +BSE results, we first notice that the transparency region is extended by 1 eV and the SHG intensity (absolute value) redistributed to higher energies. The lowest energy feature in the SHG intensity at the IPA level is a broad shoulder at about 1.3 eV. In the G 0 W 0 +BSE instead, the lowest feature in the SHG intensity is a peak at about 1.75 eV which corresponds to a two-photon resonance with the lowest optical active excitons. Overall, the SHG absolute value is lower at the G 0 W 0 +BSE than at the IPA level. The reduction is mainly due to the quasiparticle shift (a larger band gap corresponds to a lower SHG as discussed above) which is only partially compensated by the excitonic effects. Differently from the GaSe, in the InSe the excitonic effects fully compen- Arb. Units c) sate for the reduction of the SHG due to the quasiparticle shift and the local field effects. The transparency region is extended by 0.4 eV and again the SHG intensity is redistributed at higher energies and enhanced at excitonic resonances. In this case, because of the smaller quasiparticle corrections-the gap opens of 1.35 eV-the SHG at the G 0 W 0 +BSE level is slightly stronger than at the IPA level at resonance. In particular, the lowest energy feature is a broad shoulder at 1 eV within the IPA and a stronger peak at 1.32 eV corresponding to the two photons resonance of the lowest optical active exciton.

² 2 (!) ² 2 (!=2) b) G 0 W 0 + BSE Imag Real Abs 1 2 3 4 5 6 7 Energy [eV] d) ² 2 (!) ² 2 (!=2)

IV. DISCUSSION AND CONCLUSIONS

The first thing to observe is that the inclusion of manybody effects changes the SHG spectrum: both at the level of the electronic structure (quasiparticle corrections within the G 0 W 0 ) and at the level of the response function (inclusion of excitonic effects). Then, many-body effects may be important when considering technological applications (e.g. accurately determining the transparency region and excitonic resonances) and are essential when investigating the physics of excitons of these 2D systems.

On the other hand, the theoretically predicted intensities-both with and without the inclusion of many-body effects-differ substantially from the experimental values. For GaSe, Ref. 8 reports 0.7 nm/V at 0.77 eV, 1.7 nm/V at 0.92 eV and 2.4 nm/V at 1.02 eV whereas we find 0.096, 0.105, 0.111 nm/V respectively, that is a factor 7-20 smaller. For InSe, Ref. 13 reports 6.39 nm/V at 1.55 eV whereas we find 0.42 nm/V-that is a factor 15 smaller. Previously, large differences in the experimental estimate of the SHG coefficients-over a range of three orders of magnitude [START_REF] Malard | Observation of intense second harmonic generation from mos2 atomic crystals[END_REF][START_REF] Kumar | Second harmonic microscopy of monolayer mos2[END_REF][START_REF] Yin | Edge nonlinear optics on a MoS atomic monolayer[END_REF] -were observed for the MoS 2 monolayer. Eventually, Clark and coworkers [START_REF] Clark | Strong optical nonlinearity of cvd-grown mos2 monolayer as probed by wavelength-dependent secondharmonic generation[END_REF] explained that such differences depend on the model assumed to extract the value of the SHG from the experimental data. Indeed, experimentally the SHG is obtained from the measured second-harmonic intensity. Extracting the SHG coefficient is not straightforward and implies to assume a model for the system: Clark and coworkers [START_REF] Clark | Strong optical nonlinearity of cvd-grown mos2 monolayer as probed by wavelength-dependent secondharmonic generation[END_REF] argued that 2D materials should be modelled as a 'sheet', as in Refs. 4 and 12, rather than as 'bulk', as in Ref. 5. Further, they provide an expression for the difference between χ

(2) sheet and χ (2) bulk , i.e. the SHG coefficient estimated either using the 'sheet' or the 'bulk' model:

χ (2) bulk χ (2) sheet = 32π F n 2D (ω) n 2D (2ω) (n sub + 1) 3 , (2)
where n 2D and n sub are the refractive index for the 2D material and of the substrate. F is an "overall scaling factor" that depends on the transmittance and reflectance at the corresponding interfaces air-GaSe (InSe) and GaSe(InSe)-substrate as explained in the Supp. Mat. of Ref. 14. For MoS 2 , the ratio in Eq. ( 2) can be as large as 900 [START_REF] Clark | Strong optical nonlinearity of cvd-grown mos2 monolayer as probed by wavelength-dependent secondharmonic generation[END_REF] and thus accounts for the differences in the experimental estimates. [START_REF] Malard | Observation of intense second harmonic generation from mos2 atomic crystals[END_REF][START_REF] Kumar | Second harmonic microscopy of monolayer mos2[END_REF] It also indicates that the 'bulk' model leads to too large values for the SHG coefficient. This conclusion is also comforted by accurate theoretical predictions [START_REF] Grüning | Second harmonic generation in h-bn and mos2 monolayers: Role of electron-hole interaction[END_REF][START_REF] Mads | Theory of excitonic second-harmonic generation in monolayer mos2[END_REF] that reported values close to the experimental estimates that assume the 'sheet' model.

Coming back to the monochalcogenides, all the experimental values reported for the GaSe and InSe assumed a bulk model. To investigate whether this is the reason for the discrepancy with the theoretical estimates we use Eq. 2 with: n sub ≈ 1.45 (fused silica 48 ); n GaSe (ω) = n GaSe (2ω) = 2.8, 8 and n InSe (ω) = 2.86 and n InSe (2ω) = 3.35, [START_REF] Zhou | Inse monolayer: synthesis, structure and ultrahigh second-harmonic generation[END_REF] . We set F = 1, rather than the value of approximately 20 indicated in Ref. 14, since the experimental estimates already account for various gain/loss factors. We thus obtain ratios as high as ≈ 32 for GaSe and ≈ 36 for InSe. These large ratios indicate that, the experimental estimate would be significantly smaller and in substantial agreement with the theoretical predictions if using the 'sheet' model. Residual differences between the experimental and theoretical value could be due to the effects of the substrate and of the surface contributions to the SHG [START_REF] Deckoff-Jones | Observing the interplay between surface and bulk optical nonlinearities in thin van der waals crystals[END_REF] , which are neglected in the theoretical calculations. Another factor, on the side of the theoretical calculations, could be the chosen dephasing parameter as discussed in Sec. II. The so-corrected experimental estimate for GaSe is as well closer to the SHG coefficients for the bi-and trilayer, [START_REF] Jie | Layer-dependent nonlinear optical properties and stability of non-centrosymmetric modification in few-layer gase sheets[END_REF][START_REF] Tang | Layer-and frequency-dependent second harmonic generation in reflection from gase atomic crystals[END_REF] which have been obtained by modeling the few-layers as a bulk medium and accounting for interference in the multilayer system.

In summary, using accurate first-principles real-time simulations based on many-body theory, we confirm strong SHG coefficients. In particular, for InSe we obtain a SHG coefficient of ≈ 0.5nm/V at 1.67 eV-close to maximum gain and laser efficiency of Ti:Sapphire. In substantial agreement with previous theoretical works, we found that the SHG coefficients are not as strong as claimed in the experimental works of Refs. 8 and 13, differing by an order of magnitude. We argue that the experimental values are overestimated due to the assumption of a 'bulk' rather than a more appropriate 'sheet' model to extract the SHG from the measurements and they need to be reviewed. We have shown indeed that a substantial agreement between theoretical and experimental values is recovered when a 'sheet' model is used to extract the SHG from the experimental data. Note added at resubmission. Recently, we became aware of a work, Ref. 50, which reports experimental measurements of SHG in InSe in good agreement with our theoretical predictions and indeed uses the 'sheet' model to extract the SHG coefficient.

FIG. 1 .

 1 FIG. 1. [Color online] Ball-and-stick representation of the atomic structure of the GaSe (InSe) monolayer: top view (a) and lateral view (b). Gallium (Indium) atoms are in gray, Selenium in yellow.

FIG. 2 .

 2 FIG. 2. [Color online]On the left (right) the electronic band structure of InSe (GaSe) monolayer at the DFT level. Valence bands are in blue, conduction bands in green. For in-plane polarized light, the white arrow indicates the dipole-forbidden transition between the top-valence and bottom conduction bands, while the black arrow is the lowest dipole-allowed transition from the second-top valence and the bottom conduction band.

FIG. 3 .

 3 FIG. 3. [Color online] SHG [a) and b)] and optical absorption [c) and d)] of GaSe monolayer for in-plane light polarization calculated at the independent-particle [a) and c)] and at the G 0 W 0 + Bethe-Salpeter Equation level of approximation [b) and d)]. In a) and b) we plot the real part (blue line), imaginary part (red line) and the absolute value (overshadow area) of the SHG, χ[START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (ω). In c) and d) we plot the imaginary part of the dielectric function calculated at ω (continuous line) and at ω/2 (dashed line). In d) we also report the position of the two lowest bright (continuous black short vertical line) and two lowest dark (dashed black short vertical line) exciton energies for in-plane light polarization. The DFT and G 0 W 0 gap are indicated by vertical blue-violet lines.

FIG. 4 .

 4 FIG. 4. [Color online] SHG [a) and b)] and optical absorption [c) and d)] of InSe monolayer for in-plane light polarization calculated at the independent-particle [a) and c)] and at the G 0 W 0 + Bethe-Salpeter Equation level of approximation [b) and d)]. In a) and b) we plot the real part (blue line), imaginary part (red line) and the absolute value (overshadow area) of the SHG, χ[START_REF] Wang | Giant enhancement of the optical second-harmonic emission of wse2 monolayers by laser excitation at exciton resonances[END_REF] (ω). In c) and d) we plot the imaginary part of the dielectric function calculated at ω (continuous line) and at ω/2 (dashed line). In d) we also report the position of the two lowest bright (continuous black vertical line) and two lowest dark (dashed black vertical line) exciton energies for in-plane light polarization. The DFT and G 0 W 0 gap are indicated by vertical blue-violet lines

TABLE I .

 I Static limit of the SHG, χ

					G0W0+BSE)
	GaSe	110	71.0	54.0	81.2
	InSe	135	86	68.0	122

ACKNOWLEDGMENTS

The research leading to these results has received funding from the European Union Seventh Framework Program under grant agreement no. 696656 GrapheneCore1 and no. 785219 Graphene Core2. M.P. ackowledges "Tor Vergata" University for financial support through the Mission Sustainability project 2DUTOPI. E.C. acknowledges support by the Programma per Giovani Ricercatori -2014 "Rita Levi Montalcini". M.P., E.C and C.A. acknowledge PRACE for computational resources on Marconi at CINECA. C.A. acknowledges A. Zappelli for the management of the computer cluster Rosa. This publication is based upon work from COST Action TUMIEE CA17126, supported by COST (European Cooperation in Science and Technology).